0
5
0
107
78
104
82
109
80
... | ... |
@@ -221,464 +221,493 @@ |
221 | 221 |
int run(int s, Value delta = 1) { |
222 | 222 |
RangeMap<int> heap_cross_ref(_node_num, Heap::PRE_HEAP); |
223 | 223 |
Heap heap(heap_cross_ref); |
224 | 224 |
heap.push(s, 0); |
225 | 225 |
_pred[s] = -1; |
226 | 226 |
_proc_nodes.clear(); |
227 | 227 |
|
228 | 228 |
// Process nodes |
229 | 229 |
while (!heap.empty() && _excess[heap.top()] > -delta) { |
230 | 230 |
int u = heap.top(), v; |
231 | 231 |
Cost d = heap.prio() + _pi[u], dn; |
232 | 232 |
_dist[u] = heap.prio(); |
233 | 233 |
_proc_nodes.push_back(u); |
234 | 234 |
heap.pop(); |
235 | 235 |
|
236 | 236 |
// Traverse outgoing residual arcs |
237 | 237 |
int last_out = _geq ? _first_out[u+1] : _first_out[u+1] - 1; |
238 | 238 |
for (int a = _first_out[u]; a != last_out; ++a) { |
239 | 239 |
if (_res_cap[a] < delta) continue; |
240 | 240 |
v = _target[a]; |
241 | 241 |
switch (heap.state(v)) { |
242 | 242 |
case Heap::PRE_HEAP: |
243 | 243 |
heap.push(v, d + _cost[a] - _pi[v]); |
244 | 244 |
_pred[v] = a; |
245 | 245 |
break; |
246 | 246 |
case Heap::IN_HEAP: |
247 | 247 |
dn = d + _cost[a] - _pi[v]; |
248 | 248 |
if (dn < heap[v]) { |
249 | 249 |
heap.decrease(v, dn); |
250 | 250 |
_pred[v] = a; |
251 | 251 |
} |
252 | 252 |
break; |
253 | 253 |
case Heap::POST_HEAP: |
254 | 254 |
break; |
255 | 255 |
} |
256 | 256 |
} |
257 | 257 |
} |
258 | 258 |
if (heap.empty()) return -1; |
259 | 259 |
|
260 | 260 |
// Update potentials of processed nodes |
261 | 261 |
int t = heap.top(); |
262 | 262 |
Cost dt = heap.prio(); |
263 | 263 |
for (int i = 0; i < int(_proc_nodes.size()); ++i) { |
264 | 264 |
_pi[_proc_nodes[i]] += _dist[_proc_nodes[i]] - dt; |
265 | 265 |
} |
266 | 266 |
|
267 | 267 |
return t; |
268 | 268 |
} |
269 | 269 |
|
270 | 270 |
}; //class ResidualDijkstra |
271 | 271 |
|
272 | 272 |
public: |
273 | 273 |
|
274 | 274 |
/// \name Named Template Parameters |
275 | 275 |
/// @{ |
276 | 276 |
|
277 | 277 |
template <typename T> |
278 | 278 |
struct SetHeapTraits : public Traits { |
279 | 279 |
typedef T Heap; |
280 | 280 |
}; |
281 | 281 |
|
282 | 282 |
/// \brief \ref named-templ-param "Named parameter" for setting |
283 | 283 |
/// \c Heap type. |
284 | 284 |
/// |
285 | 285 |
/// \ref named-templ-param "Named parameter" for setting \c Heap |
286 | 286 |
/// type, which is used for internal Dijkstra computations. |
287 | 287 |
/// It must conform to the \ref lemon::concepts::Heap "Heap" concept, |
288 | 288 |
/// its priority type must be \c Cost and its cross reference type |
289 | 289 |
/// must be \ref RangeMap "RangeMap<int>". |
290 | 290 |
template <typename T> |
291 | 291 |
struct SetHeap |
292 | 292 |
: public CapacityScaling<GR, V, C, SetHeapTraits<T> > { |
293 | 293 |
typedef CapacityScaling<GR, V, C, SetHeapTraits<T> > Create; |
294 | 294 |
}; |
295 | 295 |
|
296 | 296 |
/// @} |
297 | 297 |
|
298 | 298 |
public: |
299 | 299 |
|
300 | 300 |
/// \brief Constructor. |
301 | 301 |
/// |
302 | 302 |
/// The constructor of the class. |
303 | 303 |
/// |
304 | 304 |
/// \param graph The digraph the algorithm runs on. |
305 | 305 |
CapacityScaling(const GR& graph) : |
306 | 306 |
_graph(graph), _node_id(graph), _arc_idf(graph), _arc_idb(graph), |
307 | 307 |
INF(std::numeric_limits<Value>::has_infinity ? |
308 | 308 |
std::numeric_limits<Value>::infinity() : |
309 | 309 |
std::numeric_limits<Value>::max()) |
310 | 310 |
{ |
311 | 311 |
// Check the number types |
312 | 312 |
LEMON_ASSERT(std::numeric_limits<Value>::is_signed, |
313 | 313 |
"The flow type of CapacityScaling must be signed"); |
314 | 314 |
LEMON_ASSERT(std::numeric_limits<Cost>::is_signed, |
315 | 315 |
"The cost type of CapacityScaling must be signed"); |
316 | 316 |
|
317 |
// Resize vectors |
|
318 |
_node_num = countNodes(_graph); |
|
319 |
_arc_num = countArcs(_graph); |
|
320 |
_res_arc_num = 2 * (_arc_num + _node_num); |
|
321 |
_root = _node_num; |
|
322 |
++_node_num; |
|
323 |
|
|
324 |
_first_out.resize(_node_num + 1); |
|
325 |
_forward.resize(_res_arc_num); |
|
326 |
_source.resize(_res_arc_num); |
|
327 |
_target.resize(_res_arc_num); |
|
328 |
_reverse.resize(_res_arc_num); |
|
329 |
|
|
330 |
_lower.resize(_res_arc_num); |
|
331 |
_upper.resize(_res_arc_num); |
|
332 |
_cost.resize(_res_arc_num); |
|
333 |
_supply.resize(_node_num); |
|
334 |
|
|
335 |
_res_cap.resize(_res_arc_num); |
|
336 |
_pi.resize(_node_num); |
|
337 |
_excess.resize(_node_num); |
|
338 |
_pred.resize(_node_num); |
|
339 |
|
|
340 |
// Copy the graph |
|
341 |
int i = 0, j = 0, k = 2 * _arc_num + _node_num - 1; |
|
342 |
for (NodeIt n(_graph); n != INVALID; ++n, ++i) { |
|
343 |
_node_id[n] = i; |
|
344 |
} |
|
345 |
i = 0; |
|
346 |
for (NodeIt n(_graph); n != INVALID; ++n, ++i) { |
|
347 |
_first_out[i] = j; |
|
348 |
for (OutArcIt a(_graph, n); a != INVALID; ++a, ++j) { |
|
349 |
_arc_idf[a] = j; |
|
350 |
_forward[j] = true; |
|
351 |
_source[j] = i; |
|
352 |
_target[j] = _node_id[_graph.runningNode(a)]; |
|
353 |
} |
|
354 |
for (InArcIt a(_graph, n); a != INVALID; ++a, ++j) { |
|
355 |
_arc_idb[a] = j; |
|
356 |
_forward[j] = false; |
|
357 |
_source[j] = i; |
|
358 |
_target[j] = _node_id[_graph.runningNode(a)]; |
|
359 |
} |
|
360 |
_forward[j] = false; |
|
361 |
_source[j] = i; |
|
362 |
_target[j] = _root; |
|
363 |
_reverse[j] = k; |
|
364 |
_forward[k] = true; |
|
365 |
_source[k] = _root; |
|
366 |
_target[k] = i; |
|
367 |
_reverse[k] = j; |
|
368 |
++j; ++k; |
|
369 |
} |
|
370 |
_first_out[i] = j; |
|
371 |
_first_out[_node_num] = k; |
|
372 |
for (ArcIt a(_graph); a != INVALID; ++a) { |
|
373 |
int fi = _arc_idf[a]; |
|
374 |
int bi = _arc_idb[a]; |
|
375 |
_reverse[fi] = bi; |
|
376 |
_reverse[bi] = fi; |
|
377 |
} |
|
378 |
|
|
379 |
// Reset |
|
317 |
// Reset data structures |
|
380 | 318 |
reset(); |
381 | 319 |
} |
382 | 320 |
|
383 | 321 |
/// \name Parameters |
384 | 322 |
/// The parameters of the algorithm can be specified using these |
385 | 323 |
/// functions. |
386 | 324 |
|
387 | 325 |
/// @{ |
388 | 326 |
|
389 | 327 |
/// \brief Set the lower bounds on the arcs. |
390 | 328 |
/// |
391 | 329 |
/// This function sets the lower bounds on the arcs. |
392 | 330 |
/// If it is not used before calling \ref run(), the lower bounds |
393 | 331 |
/// will be set to zero on all arcs. |
394 | 332 |
/// |
395 | 333 |
/// \param map An arc map storing the lower bounds. |
396 | 334 |
/// Its \c Value type must be convertible to the \c Value type |
397 | 335 |
/// of the algorithm. |
398 | 336 |
/// |
399 | 337 |
/// \return <tt>(*this)</tt> |
400 | 338 |
template <typename LowerMap> |
401 | 339 |
CapacityScaling& lowerMap(const LowerMap& map) { |
402 | 340 |
_have_lower = true; |
403 | 341 |
for (ArcIt a(_graph); a != INVALID; ++a) { |
404 | 342 |
_lower[_arc_idf[a]] = map[a]; |
405 | 343 |
_lower[_arc_idb[a]] = map[a]; |
406 | 344 |
} |
407 | 345 |
return *this; |
408 | 346 |
} |
409 | 347 |
|
410 | 348 |
/// \brief Set the upper bounds (capacities) on the arcs. |
411 | 349 |
/// |
412 | 350 |
/// This function sets the upper bounds (capacities) on the arcs. |
413 | 351 |
/// If it is not used before calling \ref run(), the upper bounds |
414 | 352 |
/// will be set to \ref INF on all arcs (i.e. the flow value will be |
415 | 353 |
/// unbounded from above). |
416 | 354 |
/// |
417 | 355 |
/// \param map An arc map storing the upper bounds. |
418 | 356 |
/// Its \c Value type must be convertible to the \c Value type |
419 | 357 |
/// of the algorithm. |
420 | 358 |
/// |
421 | 359 |
/// \return <tt>(*this)</tt> |
422 | 360 |
template<typename UpperMap> |
423 | 361 |
CapacityScaling& upperMap(const UpperMap& map) { |
424 | 362 |
for (ArcIt a(_graph); a != INVALID; ++a) { |
425 | 363 |
_upper[_arc_idf[a]] = map[a]; |
426 | 364 |
} |
427 | 365 |
return *this; |
428 | 366 |
} |
429 | 367 |
|
430 | 368 |
/// \brief Set the costs of the arcs. |
431 | 369 |
/// |
432 | 370 |
/// This function sets the costs of the arcs. |
433 | 371 |
/// If it is not used before calling \ref run(), the costs |
434 | 372 |
/// will be set to \c 1 on all arcs. |
435 | 373 |
/// |
436 | 374 |
/// \param map An arc map storing the costs. |
437 | 375 |
/// Its \c Value type must be convertible to the \c Cost type |
438 | 376 |
/// of the algorithm. |
439 | 377 |
/// |
440 | 378 |
/// \return <tt>(*this)</tt> |
441 | 379 |
template<typename CostMap> |
442 | 380 |
CapacityScaling& costMap(const CostMap& map) { |
443 | 381 |
for (ArcIt a(_graph); a != INVALID; ++a) { |
444 | 382 |
_cost[_arc_idf[a]] = map[a]; |
445 | 383 |
_cost[_arc_idb[a]] = -map[a]; |
446 | 384 |
} |
447 | 385 |
return *this; |
448 | 386 |
} |
449 | 387 |
|
450 | 388 |
/// \brief Set the supply values of the nodes. |
451 | 389 |
/// |
452 | 390 |
/// This function sets the supply values of the nodes. |
453 | 391 |
/// If neither this function nor \ref stSupply() is used before |
454 | 392 |
/// calling \ref run(), the supply of each node will be set to zero. |
455 | 393 |
/// |
456 | 394 |
/// \param map A node map storing the supply values. |
457 | 395 |
/// Its \c Value type must be convertible to the \c Value type |
458 | 396 |
/// of the algorithm. |
459 | 397 |
/// |
460 | 398 |
/// \return <tt>(*this)</tt> |
461 | 399 |
template<typename SupplyMap> |
462 | 400 |
CapacityScaling& supplyMap(const SupplyMap& map) { |
463 | 401 |
for (NodeIt n(_graph); n != INVALID; ++n) { |
464 | 402 |
_supply[_node_id[n]] = map[n]; |
465 | 403 |
} |
466 | 404 |
return *this; |
467 | 405 |
} |
468 | 406 |
|
469 | 407 |
/// \brief Set single source and target nodes and a supply value. |
470 | 408 |
/// |
471 | 409 |
/// This function sets a single source node and a single target node |
472 | 410 |
/// and the required flow value. |
473 | 411 |
/// If neither this function nor \ref supplyMap() is used before |
474 | 412 |
/// calling \ref run(), the supply of each node will be set to zero. |
475 | 413 |
/// |
476 | 414 |
/// Using this function has the same effect as using \ref supplyMap() |
477 | 415 |
/// with such a map in which \c k is assigned to \c s, \c -k is |
478 | 416 |
/// assigned to \c t and all other nodes have zero supply value. |
479 | 417 |
/// |
480 | 418 |
/// \param s The source node. |
481 | 419 |
/// \param t The target node. |
482 | 420 |
/// \param k The required amount of flow from node \c s to node \c t |
483 | 421 |
/// (i.e. the supply of \c s and the demand of \c t). |
484 | 422 |
/// |
485 | 423 |
/// \return <tt>(*this)</tt> |
486 | 424 |
CapacityScaling& stSupply(const Node& s, const Node& t, Value k) { |
487 | 425 |
for (int i = 0; i != _node_num; ++i) { |
488 | 426 |
_supply[i] = 0; |
489 | 427 |
} |
490 | 428 |
_supply[_node_id[s]] = k; |
491 | 429 |
_supply[_node_id[t]] = -k; |
492 | 430 |
return *this; |
493 | 431 |
} |
494 | 432 |
|
495 | 433 |
/// @} |
496 | 434 |
|
497 | 435 |
/// \name Execution control |
498 | 436 |
/// The algorithm can be executed using \ref run(). |
499 | 437 |
|
500 | 438 |
/// @{ |
501 | 439 |
|
502 | 440 |
/// \brief Run the algorithm. |
503 | 441 |
/// |
504 | 442 |
/// This function runs the algorithm. |
505 | 443 |
/// The paramters can be specified using functions \ref lowerMap(), |
506 | 444 |
/// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply(). |
507 | 445 |
/// For example, |
508 | 446 |
/// \code |
509 | 447 |
/// CapacityScaling<ListDigraph> cs(graph); |
510 | 448 |
/// cs.lowerMap(lower).upperMap(upper).costMap(cost) |
511 | 449 |
/// .supplyMap(sup).run(); |
512 | 450 |
/// \endcode |
513 | 451 |
/// |
514 |
/// This function can be called more than once. All the parameters |
|
515 |
/// that have been given are kept for the next call, unless |
|
516 |
/// \ref reset() is called, thus only the modified parameters |
|
517 |
/// have to be set again. See \ref reset() for examples. |
|
518 |
/// However, the underlying digraph must not be modified after this |
|
519 |
/// class have been constructed, since it copies and extends the graph. |
|
452 |
/// This function can be called more than once. All the given parameters |
|
453 |
/// are kept for the next call, unless \ref resetParams() or \ref reset() |
|
454 |
/// is used, thus only the modified parameters have to be set again. |
|
455 |
/// If the underlying digraph was also modified after the construction |
|
456 |
/// of the class (or the last \ref reset() call), then the \ref reset() |
|
457 |
/// function must be called. |
|
520 | 458 |
/// |
521 | 459 |
/// \param factor The capacity scaling factor. It must be larger than |
522 | 460 |
/// one to use scaling. If it is less or equal to one, then scaling |
523 | 461 |
/// will be disabled. |
524 | 462 |
/// |
525 | 463 |
/// \return \c INFEASIBLE if no feasible flow exists, |
526 | 464 |
/// \n \c OPTIMAL if the problem has optimal solution |
527 | 465 |
/// (i.e. it is feasible and bounded), and the algorithm has found |
528 | 466 |
/// optimal flow and node potentials (primal and dual solutions), |
529 | 467 |
/// \n \c UNBOUNDED if the digraph contains an arc of negative cost |
530 | 468 |
/// and infinite upper bound. It means that the objective function |
531 | 469 |
/// is unbounded on that arc, however, note that it could actually be |
532 | 470 |
/// bounded over the feasible flows, but this algroithm cannot handle |
533 | 471 |
/// these cases. |
534 | 472 |
/// |
535 | 473 |
/// \see ProblemType |
474 |
/// \see resetParams(), reset() |
|
536 | 475 |
ProblemType run(int factor = 4) { |
537 | 476 |
_factor = factor; |
538 | 477 |
ProblemType pt = init(); |
539 | 478 |
if (pt != OPTIMAL) return pt; |
540 | 479 |
return start(); |
541 | 480 |
} |
542 | 481 |
|
543 | 482 |
/// \brief Reset all the parameters that have been given before. |
544 | 483 |
/// |
545 | 484 |
/// This function resets all the paramaters that have been given |
546 | 485 |
/// before using functions \ref lowerMap(), \ref upperMap(), |
547 | 486 |
/// \ref costMap(), \ref supplyMap(), \ref stSupply(). |
548 | 487 |
/// |
549 |
/// It is useful for multiple run() calls. If this function is not |
|
550 |
/// used, all the parameters given before are kept for the next |
|
551 |
/// \ref run() call. |
|
552 |
/// However, the underlying digraph must not be modified after this |
|
553 |
/// |
|
488 |
/// It is useful for multiple \ref run() calls. Basically, all the given |
|
489 |
/// parameters are kept for the next \ref run() call, unless |
|
490 |
/// \ref resetParams() or \ref reset() is used. |
|
491 |
/// If the underlying digraph was also modified after the construction |
|
492 |
/// of the class or the last \ref reset() call, then the \ref reset() |
|
493 |
/// function must be used, otherwise \ref resetParams() is sufficient. |
|
554 | 494 |
/// |
555 | 495 |
/// For example, |
556 | 496 |
/// \code |
557 | 497 |
/// CapacityScaling<ListDigraph> cs(graph); |
558 | 498 |
/// |
559 | 499 |
/// // First run |
560 | 500 |
/// cs.lowerMap(lower).upperMap(upper).costMap(cost) |
561 | 501 |
/// .supplyMap(sup).run(); |
562 | 502 |
/// |
563 |
/// // Run again with modified cost map ( |
|
503 |
/// // Run again with modified cost map (resetParams() is not called, |
|
564 | 504 |
/// // so only the cost map have to be set again) |
565 | 505 |
/// cost[e] += 100; |
566 | 506 |
/// cs.costMap(cost).run(); |
567 | 507 |
/// |
568 |
/// // Run again from scratch using |
|
508 |
/// // Run again from scratch using resetParams() |
|
569 | 509 |
/// // (the lower bounds will be set to zero on all arcs) |
570 |
/// cs. |
|
510 |
/// cs.resetParams(); |
|
571 | 511 |
/// cs.upperMap(capacity).costMap(cost) |
572 | 512 |
/// .supplyMap(sup).run(); |
573 | 513 |
/// \endcode |
574 | 514 |
/// |
575 | 515 |
/// \return <tt>(*this)</tt> |
576 |
|
|
516 |
/// |
|
517 |
/// \see reset(), run() |
|
518 |
CapacityScaling& resetParams() { |
|
577 | 519 |
for (int i = 0; i != _node_num; ++i) { |
578 | 520 |
_supply[i] = 0; |
579 | 521 |
} |
580 | 522 |
for (int j = 0; j != _res_arc_num; ++j) { |
581 | 523 |
_lower[j] = 0; |
582 | 524 |
_upper[j] = INF; |
583 | 525 |
_cost[j] = _forward[j] ? 1 : -1; |
584 | 526 |
} |
585 | 527 |
_have_lower = false; |
586 | 528 |
return *this; |
587 | 529 |
} |
588 | 530 |
|
531 |
/// \brief Reset the internal data structures and all the parameters |
|
532 |
/// that have been given before. |
|
533 |
/// |
|
534 |
/// This function resets the internal data structures and all the |
|
535 |
/// paramaters that have been given before using functions \ref lowerMap(), |
|
536 |
/// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply(). |
|
537 |
/// |
|
538 |
/// It is useful for multiple \ref run() calls. Basically, all the given |
|
539 |
/// parameters are kept for the next \ref run() call, unless |
|
540 |
/// \ref resetParams() or \ref reset() is used. |
|
541 |
/// If the underlying digraph was also modified after the construction |
|
542 |
/// of the class or the last \ref reset() call, then the \ref reset() |
|
543 |
/// function must be used, otherwise \ref resetParams() is sufficient. |
|
544 |
/// |
|
545 |
/// See \ref resetParams() for examples. |
|
546 |
/// |
|
547 |
/// \return <tt>(*this)</tt> |
|
548 |
/// |
|
549 |
/// \see resetParams(), run() |
|
550 |
CapacityScaling& reset() { |
|
551 |
// Resize vectors |
|
552 |
_node_num = countNodes(_graph); |
|
553 |
_arc_num = countArcs(_graph); |
|
554 |
_res_arc_num = 2 * (_arc_num + _node_num); |
|
555 |
_root = _node_num; |
|
556 |
++_node_num; |
|
557 |
|
|
558 |
_first_out.resize(_node_num + 1); |
|
559 |
_forward.resize(_res_arc_num); |
|
560 |
_source.resize(_res_arc_num); |
|
561 |
_target.resize(_res_arc_num); |
|
562 |
_reverse.resize(_res_arc_num); |
|
563 |
|
|
564 |
_lower.resize(_res_arc_num); |
|
565 |
_upper.resize(_res_arc_num); |
|
566 |
_cost.resize(_res_arc_num); |
|
567 |
_supply.resize(_node_num); |
|
568 |
|
|
569 |
_res_cap.resize(_res_arc_num); |
|
570 |
_pi.resize(_node_num); |
|
571 |
_excess.resize(_node_num); |
|
572 |
_pred.resize(_node_num); |
|
573 |
|
|
574 |
// Copy the graph |
|
575 |
int i = 0, j = 0, k = 2 * _arc_num + _node_num - 1; |
|
576 |
for (NodeIt n(_graph); n != INVALID; ++n, ++i) { |
|
577 |
_node_id[n] = i; |
|
578 |
} |
|
579 |
i = 0; |
|
580 |
for (NodeIt n(_graph); n != INVALID; ++n, ++i) { |
|
581 |
_first_out[i] = j; |
|
582 |
for (OutArcIt a(_graph, n); a != INVALID; ++a, ++j) { |
|
583 |
_arc_idf[a] = j; |
|
584 |
_forward[j] = true; |
|
585 |
_source[j] = i; |
|
586 |
_target[j] = _node_id[_graph.runningNode(a)]; |
|
587 |
} |
|
588 |
for (InArcIt a(_graph, n); a != INVALID; ++a, ++j) { |
|
589 |
_arc_idb[a] = j; |
|
590 |
_forward[j] = false; |
|
591 |
_source[j] = i; |
|
592 |
_target[j] = _node_id[_graph.runningNode(a)]; |
|
593 |
} |
|
594 |
_forward[j] = false; |
|
595 |
_source[j] = i; |
|
596 |
_target[j] = _root; |
|
597 |
_reverse[j] = k; |
|
598 |
_forward[k] = true; |
|
599 |
_source[k] = _root; |
|
600 |
_target[k] = i; |
|
601 |
_reverse[k] = j; |
|
602 |
++j; ++k; |
|
603 |
} |
|
604 |
_first_out[i] = j; |
|
605 |
_first_out[_node_num] = k; |
|
606 |
for (ArcIt a(_graph); a != INVALID; ++a) { |
|
607 |
int fi = _arc_idf[a]; |
|
608 |
int bi = _arc_idb[a]; |
|
609 |
_reverse[fi] = bi; |
|
610 |
_reverse[bi] = fi; |
|
611 |
} |
|
612 |
|
|
613 |
// Reset parameters |
|
614 |
resetParams(); |
|
615 |
return *this; |
|
616 |
} |
|
617 |
|
|
589 | 618 |
/// @} |
590 | 619 |
|
591 | 620 |
/// \name Query Functions |
592 | 621 |
/// The results of the algorithm can be obtained using these |
593 | 622 |
/// functions.\n |
594 | 623 |
/// The \ref run() function must be called before using them. |
595 | 624 |
|
596 | 625 |
/// @{ |
597 | 626 |
|
598 | 627 |
/// \brief Return the total cost of the found flow. |
599 | 628 |
/// |
600 | 629 |
/// This function returns the total cost of the found flow. |
601 | 630 |
/// Its complexity is O(e). |
602 | 631 |
/// |
603 | 632 |
/// \note The return type of the function can be specified as a |
604 | 633 |
/// template parameter. For example, |
605 | 634 |
/// \code |
606 | 635 |
/// cs.totalCost<double>(); |
607 | 636 |
/// \endcode |
608 | 637 |
/// It is useful if the total cost cannot be stored in the \c Cost |
609 | 638 |
/// type of the algorithm, which is the default return type of the |
610 | 639 |
/// function. |
611 | 640 |
/// |
612 | 641 |
/// \pre \ref run() must be called before using this function. |
613 | 642 |
template <typename Number> |
614 | 643 |
Number totalCost() const { |
615 | 644 |
Number c = 0; |
616 | 645 |
for (ArcIt a(_graph); a != INVALID; ++a) { |
617 | 646 |
int i = _arc_idb[a]; |
618 | 647 |
c += static_cast<Number>(_res_cap[i]) * |
619 | 648 |
(-static_cast<Number>(_cost[i])); |
620 | 649 |
} |
621 | 650 |
return c; |
622 | 651 |
} |
623 | 652 |
|
624 | 653 |
#ifndef DOXYGEN |
625 | 654 |
Cost totalCost() const { |
626 | 655 |
return totalCost<Cost>(); |
627 | 656 |
} |
628 | 657 |
#endif |
629 | 658 |
|
630 | 659 |
/// \brief Return the flow on the given arc. |
631 | 660 |
/// |
632 | 661 |
/// This function returns the flow on the given arc. |
633 | 662 |
/// |
634 | 663 |
/// \pre \ref run() must be called before using this function. |
635 | 664 |
Value flow(const Arc& a) const { |
636 | 665 |
return _res_cap[_arc_idb[a]]; |
637 | 666 |
} |
638 | 667 |
|
639 | 668 |
/// \brief Return the flow map (the primal solution). |
640 | 669 |
/// |
641 | 670 |
/// This function copies the flow value on each arc into the given |
642 | 671 |
/// map. The \c Value type of the algorithm must be convertible to |
643 | 672 |
/// the \c Value type of the map. |
644 | 673 |
/// |
645 | 674 |
/// \pre \ref run() must be called before using this function. |
646 | 675 |
template <typename FlowMap> |
647 | 676 |
void flowMap(FlowMap &map) const { |
648 | 677 |
for (ArcIt a(_graph); a != INVALID; ++a) { |
649 | 678 |
map.set(a, _res_cap[_arc_idb[a]]); |
650 | 679 |
} |
651 | 680 |
} |
652 | 681 |
|
653 | 682 |
/// \brief Return the potential (dual value) of the given node. |
654 | 683 |
/// |
655 | 684 |
/// This function returns the potential (dual value) of the |
656 | 685 |
/// given node. |
657 | 686 |
/// |
658 | 687 |
/// \pre \ref run() must be called before using this function. |
659 | 688 |
Cost potential(const Node& n) const { |
660 | 689 |
return _pi[_node_id[n]]; |
661 | 690 |
} |
662 | 691 |
|
663 | 692 |
/// \brief Return the potential map (the dual solution). |
664 | 693 |
/// |
665 | 694 |
/// This function copies the potential (dual value) of each node |
666 | 695 |
/// into the given map. |
667 | 696 |
/// The \c Cost type of the algorithm must be convertible to the |
668 | 697 |
/// \c Value type of the map. |
669 | 698 |
/// |
670 | 699 |
/// \pre \ref run() must be called before using this function. |
671 | 700 |
template <typename PotentialMap> |
672 | 701 |
void potentialMap(PotentialMap &map) const { |
673 | 702 |
for (NodeIt n(_graph); n != INVALID; ++n) { |
674 | 703 |
map.set(n, _pi[_node_id[n]]); |
675 | 704 |
} |
676 | 705 |
} |
677 | 706 |
|
678 | 707 |
/// @} |
679 | 708 |
|
680 | 709 |
private: |
681 | 710 |
|
682 | 711 |
// Initialize the algorithm |
683 | 712 |
ProblemType init() { |
684 | 713 |
if (_node_num <= 1) return INFEASIBLE; |
... | ... |
@@ -240,476 +240,498 @@ |
240 | 240 |
int _res_node_num; |
241 | 241 |
int _res_arc_num; |
242 | 242 |
int _root; |
243 | 243 |
|
244 | 244 |
// Parameters of the problem |
245 | 245 |
bool _have_lower; |
246 | 246 |
Value _sum_supply; |
247 | 247 |
|
248 | 248 |
// Data structures for storing the digraph |
249 | 249 |
IntNodeMap _node_id; |
250 | 250 |
IntArcMap _arc_idf; |
251 | 251 |
IntArcMap _arc_idb; |
252 | 252 |
IntVector _first_out; |
253 | 253 |
BoolVector _forward; |
254 | 254 |
IntVector _source; |
255 | 255 |
IntVector _target; |
256 | 256 |
IntVector _reverse; |
257 | 257 |
|
258 | 258 |
// Node and arc data |
259 | 259 |
ValueVector _lower; |
260 | 260 |
ValueVector _upper; |
261 | 261 |
CostVector _scost; |
262 | 262 |
ValueVector _supply; |
263 | 263 |
|
264 | 264 |
ValueVector _res_cap; |
265 | 265 |
LargeCostVector _cost; |
266 | 266 |
LargeCostVector _pi; |
267 | 267 |
ValueVector _excess; |
268 | 268 |
IntVector _next_out; |
269 | 269 |
std::deque<int> _active_nodes; |
270 | 270 |
|
271 | 271 |
// Data for scaling |
272 | 272 |
LargeCost _epsilon; |
273 | 273 |
int _alpha; |
274 | 274 |
|
275 | 275 |
// Data for a StaticDigraph structure |
276 | 276 |
typedef std::pair<int, int> IntPair; |
277 | 277 |
StaticDigraph _sgr; |
278 | 278 |
std::vector<IntPair> _arc_vec; |
279 | 279 |
std::vector<LargeCost> _cost_vec; |
280 | 280 |
LargeCostArcMap _cost_map; |
281 | 281 |
LargeCostNodeMap _pi_map; |
282 | 282 |
|
283 | 283 |
public: |
284 | 284 |
|
285 | 285 |
/// \brief Constant for infinite upper bounds (capacities). |
286 | 286 |
/// |
287 | 287 |
/// Constant for infinite upper bounds (capacities). |
288 | 288 |
/// It is \c std::numeric_limits<Value>::infinity() if available, |
289 | 289 |
/// \c std::numeric_limits<Value>::max() otherwise. |
290 | 290 |
const Value INF; |
291 | 291 |
|
292 | 292 |
public: |
293 | 293 |
|
294 | 294 |
/// \name Named Template Parameters |
295 | 295 |
/// @{ |
296 | 296 |
|
297 | 297 |
template <typename T> |
298 | 298 |
struct SetLargeCostTraits : public Traits { |
299 | 299 |
typedef T LargeCost; |
300 | 300 |
}; |
301 | 301 |
|
302 | 302 |
/// \brief \ref named-templ-param "Named parameter" for setting |
303 | 303 |
/// \c LargeCost type. |
304 | 304 |
/// |
305 | 305 |
/// \ref named-templ-param "Named parameter" for setting \c LargeCost |
306 | 306 |
/// type, which is used for internal computations in the algorithm. |
307 | 307 |
/// \c Cost must be convertible to \c LargeCost. |
308 | 308 |
template <typename T> |
309 | 309 |
struct SetLargeCost |
310 | 310 |
: public CostScaling<GR, V, C, SetLargeCostTraits<T> > { |
311 | 311 |
typedef CostScaling<GR, V, C, SetLargeCostTraits<T> > Create; |
312 | 312 |
}; |
313 | 313 |
|
314 | 314 |
/// @} |
315 | 315 |
|
316 | 316 |
public: |
317 | 317 |
|
318 | 318 |
/// \brief Constructor. |
319 | 319 |
/// |
320 | 320 |
/// The constructor of the class. |
321 | 321 |
/// |
322 | 322 |
/// \param graph The digraph the algorithm runs on. |
323 | 323 |
CostScaling(const GR& graph) : |
324 | 324 |
_graph(graph), _node_id(graph), _arc_idf(graph), _arc_idb(graph), |
325 | 325 |
_cost_map(_cost_vec), _pi_map(_pi), |
326 | 326 |
INF(std::numeric_limits<Value>::has_infinity ? |
327 | 327 |
std::numeric_limits<Value>::infinity() : |
328 | 328 |
std::numeric_limits<Value>::max()) |
329 | 329 |
{ |
330 | 330 |
// Check the number types |
331 | 331 |
LEMON_ASSERT(std::numeric_limits<Value>::is_signed, |
332 | 332 |
"The flow type of CostScaling must be signed"); |
333 | 333 |
LEMON_ASSERT(std::numeric_limits<Cost>::is_signed, |
334 | 334 |
"The cost type of CostScaling must be signed"); |
335 | 335 |
|
336 |
// Resize vectors |
|
337 |
_node_num = countNodes(_graph); |
|
338 |
_arc_num = countArcs(_graph); |
|
339 |
_res_node_num = _node_num + 1; |
|
340 |
_res_arc_num = 2 * (_arc_num + _node_num); |
|
341 |
_root = _node_num; |
|
342 |
|
|
343 |
_first_out.resize(_res_node_num + 1); |
|
344 |
_forward.resize(_res_arc_num); |
|
345 |
_source.resize(_res_arc_num); |
|
346 |
_target.resize(_res_arc_num); |
|
347 |
_reverse.resize(_res_arc_num); |
|
348 |
|
|
349 |
_lower.resize(_res_arc_num); |
|
350 |
_upper.resize(_res_arc_num); |
|
351 |
_scost.resize(_res_arc_num); |
|
352 |
_supply.resize(_res_node_num); |
|
353 |
|
|
354 |
_res_cap.resize(_res_arc_num); |
|
355 |
_cost.resize(_res_arc_num); |
|
356 |
_pi.resize(_res_node_num); |
|
357 |
_excess.resize(_res_node_num); |
|
358 |
_next_out.resize(_res_node_num); |
|
359 |
|
|
360 |
_arc_vec.reserve(_res_arc_num); |
|
361 |
_cost_vec.reserve(_res_arc_num); |
|
362 |
|
|
363 |
// Copy the graph |
|
364 |
int i = 0, j = 0, k = 2 * _arc_num + _node_num; |
|
365 |
for (NodeIt n(_graph); n != INVALID; ++n, ++i) { |
|
366 |
_node_id[n] = i; |
|
367 |
} |
|
368 |
i = 0; |
|
369 |
for (NodeIt n(_graph); n != INVALID; ++n, ++i) { |
|
370 |
_first_out[i] = j; |
|
371 |
for (OutArcIt a(_graph, n); a != INVALID; ++a, ++j) { |
|
372 |
_arc_idf[a] = j; |
|
373 |
_forward[j] = true; |
|
374 |
_source[j] = i; |
|
375 |
_target[j] = _node_id[_graph.runningNode(a)]; |
|
376 |
} |
|
377 |
for (InArcIt a(_graph, n); a != INVALID; ++a, ++j) { |
|
378 |
_arc_idb[a] = j; |
|
379 |
_forward[j] = false; |
|
380 |
_source[j] = i; |
|
381 |
_target[j] = _node_id[_graph.runningNode(a)]; |
|
382 |
} |
|
383 |
_forward[j] = false; |
|
384 |
_source[j] = i; |
|
385 |
_target[j] = _root; |
|
386 |
_reverse[j] = k; |
|
387 |
_forward[k] = true; |
|
388 |
_source[k] = _root; |
|
389 |
_target[k] = i; |
|
390 |
_reverse[k] = j; |
|
391 |
++j; ++k; |
|
392 |
} |
|
393 |
_first_out[i] = j; |
|
394 |
_first_out[_res_node_num] = k; |
|
395 |
for (ArcIt a(_graph); a != INVALID; ++a) { |
|
396 |
int fi = _arc_idf[a]; |
|
397 |
int bi = _arc_idb[a]; |
|
398 |
_reverse[fi] = bi; |
|
399 |
_reverse[bi] = fi; |
|
400 |
} |
|
401 |
|
|
402 |
// Reset |
|
336 |
// Reset data structures |
|
403 | 337 |
reset(); |
404 | 338 |
} |
405 | 339 |
|
406 | 340 |
/// \name Parameters |
407 | 341 |
/// The parameters of the algorithm can be specified using these |
408 | 342 |
/// functions. |
409 | 343 |
|
410 | 344 |
/// @{ |
411 | 345 |
|
412 | 346 |
/// \brief Set the lower bounds on the arcs. |
413 | 347 |
/// |
414 | 348 |
/// This function sets the lower bounds on the arcs. |
415 | 349 |
/// If it is not used before calling \ref run(), the lower bounds |
416 | 350 |
/// will be set to zero on all arcs. |
417 | 351 |
/// |
418 | 352 |
/// \param map An arc map storing the lower bounds. |
419 | 353 |
/// Its \c Value type must be convertible to the \c Value type |
420 | 354 |
/// of the algorithm. |
421 | 355 |
/// |
422 | 356 |
/// \return <tt>(*this)</tt> |
423 | 357 |
template <typename LowerMap> |
424 | 358 |
CostScaling& lowerMap(const LowerMap& map) { |
425 | 359 |
_have_lower = true; |
426 | 360 |
for (ArcIt a(_graph); a != INVALID; ++a) { |
427 | 361 |
_lower[_arc_idf[a]] = map[a]; |
428 | 362 |
_lower[_arc_idb[a]] = map[a]; |
429 | 363 |
} |
430 | 364 |
return *this; |
431 | 365 |
} |
432 | 366 |
|
433 | 367 |
/// \brief Set the upper bounds (capacities) on the arcs. |
434 | 368 |
/// |
435 | 369 |
/// This function sets the upper bounds (capacities) on the arcs. |
436 | 370 |
/// If it is not used before calling \ref run(), the upper bounds |
437 | 371 |
/// will be set to \ref INF on all arcs (i.e. the flow value will be |
438 | 372 |
/// unbounded from above). |
439 | 373 |
/// |
440 | 374 |
/// \param map An arc map storing the upper bounds. |
441 | 375 |
/// Its \c Value type must be convertible to the \c Value type |
442 | 376 |
/// of the algorithm. |
443 | 377 |
/// |
444 | 378 |
/// \return <tt>(*this)</tt> |
445 | 379 |
template<typename UpperMap> |
446 | 380 |
CostScaling& upperMap(const UpperMap& map) { |
447 | 381 |
for (ArcIt a(_graph); a != INVALID; ++a) { |
448 | 382 |
_upper[_arc_idf[a]] = map[a]; |
449 | 383 |
} |
450 | 384 |
return *this; |
451 | 385 |
} |
452 | 386 |
|
453 | 387 |
/// \brief Set the costs of the arcs. |
454 | 388 |
/// |
455 | 389 |
/// This function sets the costs of the arcs. |
456 | 390 |
/// If it is not used before calling \ref run(), the costs |
457 | 391 |
/// will be set to \c 1 on all arcs. |
458 | 392 |
/// |
459 | 393 |
/// \param map An arc map storing the costs. |
460 | 394 |
/// Its \c Value type must be convertible to the \c Cost type |
461 | 395 |
/// of the algorithm. |
462 | 396 |
/// |
463 | 397 |
/// \return <tt>(*this)</tt> |
464 | 398 |
template<typename CostMap> |
465 | 399 |
CostScaling& costMap(const CostMap& map) { |
466 | 400 |
for (ArcIt a(_graph); a != INVALID; ++a) { |
467 | 401 |
_scost[_arc_idf[a]] = map[a]; |
468 | 402 |
_scost[_arc_idb[a]] = -map[a]; |
469 | 403 |
} |
470 | 404 |
return *this; |
471 | 405 |
} |
472 | 406 |
|
473 | 407 |
/// \brief Set the supply values of the nodes. |
474 | 408 |
/// |
475 | 409 |
/// This function sets the supply values of the nodes. |
476 | 410 |
/// If neither this function nor \ref stSupply() is used before |
477 | 411 |
/// calling \ref run(), the supply of each node will be set to zero. |
478 | 412 |
/// |
479 | 413 |
/// \param map A node map storing the supply values. |
480 | 414 |
/// Its \c Value type must be convertible to the \c Value type |
481 | 415 |
/// of the algorithm. |
482 | 416 |
/// |
483 | 417 |
/// \return <tt>(*this)</tt> |
484 | 418 |
template<typename SupplyMap> |
485 | 419 |
CostScaling& supplyMap(const SupplyMap& map) { |
486 | 420 |
for (NodeIt n(_graph); n != INVALID; ++n) { |
487 | 421 |
_supply[_node_id[n]] = map[n]; |
488 | 422 |
} |
489 | 423 |
return *this; |
490 | 424 |
} |
491 | 425 |
|
492 | 426 |
/// \brief Set single source and target nodes and a supply value. |
493 | 427 |
/// |
494 | 428 |
/// This function sets a single source node and a single target node |
495 | 429 |
/// and the required flow value. |
496 | 430 |
/// If neither this function nor \ref supplyMap() is used before |
497 | 431 |
/// calling \ref run(), the supply of each node will be set to zero. |
498 | 432 |
/// |
499 | 433 |
/// Using this function has the same effect as using \ref supplyMap() |
500 | 434 |
/// with such a map in which \c k is assigned to \c s, \c -k is |
501 | 435 |
/// assigned to \c t and all other nodes have zero supply value. |
502 | 436 |
/// |
503 | 437 |
/// \param s The source node. |
504 | 438 |
/// \param t The target node. |
505 | 439 |
/// \param k The required amount of flow from node \c s to node \c t |
506 | 440 |
/// (i.e. the supply of \c s and the demand of \c t). |
507 | 441 |
/// |
508 | 442 |
/// \return <tt>(*this)</tt> |
509 | 443 |
CostScaling& stSupply(const Node& s, const Node& t, Value k) { |
510 | 444 |
for (int i = 0; i != _res_node_num; ++i) { |
511 | 445 |
_supply[i] = 0; |
512 | 446 |
} |
513 | 447 |
_supply[_node_id[s]] = k; |
514 | 448 |
_supply[_node_id[t]] = -k; |
515 | 449 |
return *this; |
516 | 450 |
} |
517 | 451 |
|
518 | 452 |
/// @} |
519 | 453 |
|
520 | 454 |
/// \name Execution control |
521 | 455 |
/// The algorithm can be executed using \ref run(). |
522 | 456 |
|
523 | 457 |
/// @{ |
524 | 458 |
|
525 | 459 |
/// \brief Run the algorithm. |
526 | 460 |
/// |
527 | 461 |
/// This function runs the algorithm. |
528 | 462 |
/// The paramters can be specified using functions \ref lowerMap(), |
529 | 463 |
/// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply(). |
530 | 464 |
/// For example, |
531 | 465 |
/// \code |
532 | 466 |
/// CostScaling<ListDigraph> cs(graph); |
533 | 467 |
/// cs.lowerMap(lower).upperMap(upper).costMap(cost) |
534 | 468 |
/// .supplyMap(sup).run(); |
535 | 469 |
/// \endcode |
536 | 470 |
/// |
537 |
/// This function can be called more than once. All the parameters |
|
538 |
/// that have been given are kept for the next call, unless |
|
539 |
/// \ref reset() is called, thus only the modified parameters |
|
540 |
/// have to be set again. See \ref reset() for examples. |
|
541 |
/// However, the underlying digraph must not be modified after this |
|
542 |
/// class have been constructed, since it copies and extends the graph. |
|
471 |
/// This function can be called more than once. All the given parameters |
|
472 |
/// are kept for the next call, unless \ref resetParams() or \ref reset() |
|
473 |
/// is used, thus only the modified parameters have to be set again. |
|
474 |
/// If the underlying digraph was also modified after the construction |
|
475 |
/// of the class (or the last \ref reset() call), then the \ref reset() |
|
476 |
/// function must be called. |
|
543 | 477 |
/// |
544 | 478 |
/// \param method The internal method that will be used in the |
545 | 479 |
/// algorithm. For more information, see \ref Method. |
546 | 480 |
/// \param factor The cost scaling factor. It must be larger than one. |
547 | 481 |
/// |
548 | 482 |
/// \return \c INFEASIBLE if no feasible flow exists, |
549 | 483 |
/// \n \c OPTIMAL if the problem has optimal solution |
550 | 484 |
/// (i.e. it is feasible and bounded), and the algorithm has found |
551 | 485 |
/// optimal flow and node potentials (primal and dual solutions), |
552 | 486 |
/// \n \c UNBOUNDED if the digraph contains an arc of negative cost |
553 | 487 |
/// and infinite upper bound. It means that the objective function |
554 | 488 |
/// is unbounded on that arc, however, note that it could actually be |
555 | 489 |
/// bounded over the feasible flows, but this algroithm cannot handle |
556 | 490 |
/// these cases. |
557 | 491 |
/// |
558 | 492 |
/// \see ProblemType, Method |
493 |
/// \see resetParams(), reset() |
|
559 | 494 |
ProblemType run(Method method = PARTIAL_AUGMENT, int factor = 8) { |
560 | 495 |
_alpha = factor; |
561 | 496 |
ProblemType pt = init(); |
562 | 497 |
if (pt != OPTIMAL) return pt; |
563 | 498 |
start(method); |
564 | 499 |
return OPTIMAL; |
565 | 500 |
} |
566 | 501 |
|
567 | 502 |
/// \brief Reset all the parameters that have been given before. |
568 | 503 |
/// |
569 | 504 |
/// This function resets all the paramaters that have been given |
570 | 505 |
/// before using functions \ref lowerMap(), \ref upperMap(), |
571 | 506 |
/// \ref costMap(), \ref supplyMap(), \ref stSupply(). |
572 | 507 |
/// |
573 |
/// It is useful for multiple run() calls. If this function is not |
|
574 |
/// used, all the parameters given before are kept for the next |
|
575 |
/// \ref run() call. |
|
576 |
/// However, the underlying digraph must not be modified after this |
|
577 |
/// |
|
508 |
/// It is useful for multiple \ref run() calls. Basically, all the given |
|
509 |
/// parameters are kept for the next \ref run() call, unless |
|
510 |
/// \ref resetParams() or \ref reset() is used. |
|
511 |
/// If the underlying digraph was also modified after the construction |
|
512 |
/// of the class or the last \ref reset() call, then the \ref reset() |
|
513 |
/// function must be used, otherwise \ref resetParams() is sufficient. |
|
578 | 514 |
/// |
579 | 515 |
/// For example, |
580 | 516 |
/// \code |
581 | 517 |
/// CostScaling<ListDigraph> cs(graph); |
582 | 518 |
/// |
583 | 519 |
/// // First run |
584 | 520 |
/// cs.lowerMap(lower).upperMap(upper).costMap(cost) |
585 | 521 |
/// .supplyMap(sup).run(); |
586 | 522 |
/// |
587 |
/// // Run again with modified cost map ( |
|
523 |
/// // Run again with modified cost map (resetParams() is not called, |
|
588 | 524 |
/// // so only the cost map have to be set again) |
589 | 525 |
/// cost[e] += 100; |
590 | 526 |
/// cs.costMap(cost).run(); |
591 | 527 |
/// |
592 |
/// // Run again from scratch using |
|
528 |
/// // Run again from scratch using resetParams() |
|
593 | 529 |
/// // (the lower bounds will be set to zero on all arcs) |
594 |
/// cs. |
|
530 |
/// cs.resetParams(); |
|
595 | 531 |
/// cs.upperMap(capacity).costMap(cost) |
596 | 532 |
/// .supplyMap(sup).run(); |
597 | 533 |
/// \endcode |
598 | 534 |
/// |
599 | 535 |
/// \return <tt>(*this)</tt> |
600 |
|
|
536 |
/// |
|
537 |
/// \see reset(), run() |
|
538 |
CostScaling& resetParams() { |
|
601 | 539 |
for (int i = 0; i != _res_node_num; ++i) { |
602 | 540 |
_supply[i] = 0; |
603 | 541 |
} |
604 | 542 |
int limit = _first_out[_root]; |
605 | 543 |
for (int j = 0; j != limit; ++j) { |
606 | 544 |
_lower[j] = 0; |
607 | 545 |
_upper[j] = INF; |
608 | 546 |
_scost[j] = _forward[j] ? 1 : -1; |
609 | 547 |
} |
610 | 548 |
for (int j = limit; j != _res_arc_num; ++j) { |
611 | 549 |
_lower[j] = 0; |
612 | 550 |
_upper[j] = INF; |
613 | 551 |
_scost[j] = 0; |
614 | 552 |
_scost[_reverse[j]] = 0; |
615 | 553 |
} |
616 | 554 |
_have_lower = false; |
617 | 555 |
return *this; |
618 | 556 |
} |
619 | 557 |
|
558 |
/// \brief Reset all the parameters that have been given before. |
|
559 |
/// |
|
560 |
/// This function resets all the paramaters that have been given |
|
561 |
/// before using functions \ref lowerMap(), \ref upperMap(), |
|
562 |
/// \ref costMap(), \ref supplyMap(), \ref stSupply(). |
|
563 |
/// |
|
564 |
/// It is useful for multiple run() calls. If this function is not |
|
565 |
/// used, all the parameters given before are kept for the next |
|
566 |
/// \ref run() call. |
|
567 |
/// However, the underlying digraph must not be modified after this |
|
568 |
/// class have been constructed, since it copies and extends the graph. |
|
569 |
/// \return <tt>(*this)</tt> |
|
570 |
CostScaling& reset() { |
|
571 |
// Resize vectors |
|
572 |
_node_num = countNodes(_graph); |
|
573 |
_arc_num = countArcs(_graph); |
|
574 |
_res_node_num = _node_num + 1; |
|
575 |
_res_arc_num = 2 * (_arc_num + _node_num); |
|
576 |
_root = _node_num; |
|
577 |
|
|
578 |
_first_out.resize(_res_node_num + 1); |
|
579 |
_forward.resize(_res_arc_num); |
|
580 |
_source.resize(_res_arc_num); |
|
581 |
_target.resize(_res_arc_num); |
|
582 |
_reverse.resize(_res_arc_num); |
|
583 |
|
|
584 |
_lower.resize(_res_arc_num); |
|
585 |
_upper.resize(_res_arc_num); |
|
586 |
_scost.resize(_res_arc_num); |
|
587 |
_supply.resize(_res_node_num); |
|
588 |
|
|
589 |
_res_cap.resize(_res_arc_num); |
|
590 |
_cost.resize(_res_arc_num); |
|
591 |
_pi.resize(_res_node_num); |
|
592 |
_excess.resize(_res_node_num); |
|
593 |
_next_out.resize(_res_node_num); |
|
594 |
|
|
595 |
_arc_vec.reserve(_res_arc_num); |
|
596 |
_cost_vec.reserve(_res_arc_num); |
|
597 |
|
|
598 |
// Copy the graph |
|
599 |
int i = 0, j = 0, k = 2 * _arc_num + _node_num; |
|
600 |
for (NodeIt n(_graph); n != INVALID; ++n, ++i) { |
|
601 |
_node_id[n] = i; |
|
602 |
} |
|
603 |
i = 0; |
|
604 |
for (NodeIt n(_graph); n != INVALID; ++n, ++i) { |
|
605 |
_first_out[i] = j; |
|
606 |
for (OutArcIt a(_graph, n); a != INVALID; ++a, ++j) { |
|
607 |
_arc_idf[a] = j; |
|
608 |
_forward[j] = true; |
|
609 |
_source[j] = i; |
|
610 |
_target[j] = _node_id[_graph.runningNode(a)]; |
|
611 |
} |
|
612 |
for (InArcIt a(_graph, n); a != INVALID; ++a, ++j) { |
|
613 |
_arc_idb[a] = j; |
|
614 |
_forward[j] = false; |
|
615 |
_source[j] = i; |
|
616 |
_target[j] = _node_id[_graph.runningNode(a)]; |
|
617 |
} |
|
618 |
_forward[j] = false; |
|
619 |
_source[j] = i; |
|
620 |
_target[j] = _root; |
|
621 |
_reverse[j] = k; |
|
622 |
_forward[k] = true; |
|
623 |
_source[k] = _root; |
|
624 |
_target[k] = i; |
|
625 |
_reverse[k] = j; |
|
626 |
++j; ++k; |
|
627 |
} |
|
628 |
_first_out[i] = j; |
|
629 |
_first_out[_res_node_num] = k; |
|
630 |
for (ArcIt a(_graph); a != INVALID; ++a) { |
|
631 |
int fi = _arc_idf[a]; |
|
632 |
int bi = _arc_idb[a]; |
|
633 |
_reverse[fi] = bi; |
|
634 |
_reverse[bi] = fi; |
|
635 |
} |
|
636 |
|
|
637 |
// Reset parameters |
|
638 |
resetParams(); |
|
639 |
return *this; |
|
640 |
} |
|
641 |
|
|
620 | 642 |
/// @} |
621 | 643 |
|
622 | 644 |
/// \name Query Functions |
623 | 645 |
/// The results of the algorithm can be obtained using these |
624 | 646 |
/// functions.\n |
625 | 647 |
/// The \ref run() function must be called before using them. |
626 | 648 |
|
627 | 649 |
/// @{ |
628 | 650 |
|
629 | 651 |
/// \brief Return the total cost of the found flow. |
630 | 652 |
/// |
631 | 653 |
/// This function returns the total cost of the found flow. |
632 | 654 |
/// Its complexity is O(e). |
633 | 655 |
/// |
634 | 656 |
/// \note The return type of the function can be specified as a |
635 | 657 |
/// template parameter. For example, |
636 | 658 |
/// \code |
637 | 659 |
/// cs.totalCost<double>(); |
638 | 660 |
/// \endcode |
639 | 661 |
/// It is useful if the total cost cannot be stored in the \c Cost |
640 | 662 |
/// type of the algorithm, which is the default return type of the |
641 | 663 |
/// function. |
642 | 664 |
/// |
643 | 665 |
/// \pre \ref run() must be called before using this function. |
644 | 666 |
template <typename Number> |
645 | 667 |
Number totalCost() const { |
646 | 668 |
Number c = 0; |
647 | 669 |
for (ArcIt a(_graph); a != INVALID; ++a) { |
648 | 670 |
int i = _arc_idb[a]; |
649 | 671 |
c += static_cast<Number>(_res_cap[i]) * |
650 | 672 |
(-static_cast<Number>(_scost[i])); |
651 | 673 |
} |
652 | 674 |
return c; |
653 | 675 |
} |
654 | 676 |
|
655 | 677 |
#ifndef DOXYGEN |
656 | 678 |
Cost totalCost() const { |
657 | 679 |
return totalCost<Cost>(); |
658 | 680 |
} |
659 | 681 |
#endif |
660 | 682 |
|
661 | 683 |
/// \brief Return the flow on the given arc. |
662 | 684 |
/// |
663 | 685 |
/// This function returns the flow on the given arc. |
664 | 686 |
/// |
665 | 687 |
/// \pre \ref run() must be called before using this function. |
666 | 688 |
Value flow(const Arc& a) const { |
667 | 689 |
return _res_cap[_arc_idb[a]]; |
668 | 690 |
} |
669 | 691 |
|
670 | 692 |
/// \brief Return the flow map (the primal solution). |
671 | 693 |
/// |
672 | 694 |
/// This function copies the flow value on each arc into the given |
673 | 695 |
/// map. The \c Value type of the algorithm must be convertible to |
674 | 696 |
/// the \c Value type of the map. |
675 | 697 |
/// |
676 | 698 |
/// \pre \ref run() must be called before using this function. |
677 | 699 |
template <typename FlowMap> |
678 | 700 |
void flowMap(FlowMap &map) const { |
679 | 701 |
for (ArcIt a(_graph); a != INVALID; ++a) { |
680 | 702 |
map.set(a, _res_cap[_arc_idb[a]]); |
681 | 703 |
} |
682 | 704 |
} |
683 | 705 |
|
684 | 706 |
/// \brief Return the potential (dual value) of the given node. |
685 | 707 |
/// |
686 | 708 |
/// This function returns the potential (dual value) of the |
687 | 709 |
/// given node. |
688 | 710 |
/// |
689 | 711 |
/// \pre \ref run() must be called before using this function. |
690 | 712 |
Cost potential(const Node& n) const { |
691 | 713 |
return static_cast<Cost>(_pi[_node_id[n]]); |
692 | 714 |
} |
693 | 715 |
|
694 | 716 |
/// \brief Return the potential map (the dual solution). |
695 | 717 |
/// |
696 | 718 |
/// This function copies the potential (dual value) of each node |
697 | 719 |
/// into the given map. |
698 | 720 |
/// The \c Cost type of the algorithm must be convertible to the |
699 | 721 |
/// \c Value type of the map. |
700 | 722 |
/// |
701 | 723 |
/// \pre \ref run() must be called before using this function. |
702 | 724 |
template <typename PotentialMap> |
703 | 725 |
void potentialMap(PotentialMap &map) const { |
704 | 726 |
for (NodeIt n(_graph); n != INVALID; ++n) { |
705 | 727 |
map.set(n, static_cast<Cost>(_pi[_node_id[n]])); |
706 | 728 |
} |
707 | 729 |
} |
708 | 730 |
|
709 | 731 |
/// @} |
710 | 732 |
|
711 | 733 |
private: |
712 | 734 |
|
713 | 735 |
// Initialize the algorithm |
714 | 736 |
ProblemType init() { |
715 | 737 |
if (_res_node_num <= 1) return INFEASIBLE; |
... | ... |
@@ -157,472 +157,501 @@ |
157 | 157 |
typedef KT Key; |
158 | 158 |
typedef VT Value; |
159 | 159 |
|
160 | 160 |
StaticVectorMap(std::vector<Value>& v) : _v(v) {} |
161 | 161 |
|
162 | 162 |
const Value& operator[](const Key& key) const { |
163 | 163 |
return _v[StaticDigraph::id(key)]; |
164 | 164 |
} |
165 | 165 |
|
166 | 166 |
Value& operator[](const Key& key) { |
167 | 167 |
return _v[StaticDigraph::id(key)]; |
168 | 168 |
} |
169 | 169 |
|
170 | 170 |
void set(const Key& key, const Value& val) { |
171 | 171 |
_v[StaticDigraph::id(key)] = val; |
172 | 172 |
} |
173 | 173 |
|
174 | 174 |
private: |
175 | 175 |
std::vector<Value>& _v; |
176 | 176 |
}; |
177 | 177 |
|
178 | 178 |
typedef StaticVectorMap<StaticDigraph::Node, Cost> CostNodeMap; |
179 | 179 |
typedef StaticVectorMap<StaticDigraph::Arc, Cost> CostArcMap; |
180 | 180 |
|
181 | 181 |
private: |
182 | 182 |
|
183 | 183 |
|
184 | 184 |
// Data related to the underlying digraph |
185 | 185 |
const GR &_graph; |
186 | 186 |
int _node_num; |
187 | 187 |
int _arc_num; |
188 | 188 |
int _res_node_num; |
189 | 189 |
int _res_arc_num; |
190 | 190 |
int _root; |
191 | 191 |
|
192 | 192 |
// Parameters of the problem |
193 | 193 |
bool _have_lower; |
194 | 194 |
Value _sum_supply; |
195 | 195 |
|
196 | 196 |
// Data structures for storing the digraph |
197 | 197 |
IntNodeMap _node_id; |
198 | 198 |
IntArcMap _arc_idf; |
199 | 199 |
IntArcMap _arc_idb; |
200 | 200 |
IntVector _first_out; |
201 | 201 |
CharVector _forward; |
202 | 202 |
IntVector _source; |
203 | 203 |
IntVector _target; |
204 | 204 |
IntVector _reverse; |
205 | 205 |
|
206 | 206 |
// Node and arc data |
207 | 207 |
ValueVector _lower; |
208 | 208 |
ValueVector _upper; |
209 | 209 |
CostVector _cost; |
210 | 210 |
ValueVector _supply; |
211 | 211 |
|
212 | 212 |
ValueVector _res_cap; |
213 | 213 |
CostVector _pi; |
214 | 214 |
|
215 | 215 |
// Data for a StaticDigraph structure |
216 | 216 |
typedef std::pair<int, int> IntPair; |
217 | 217 |
StaticDigraph _sgr; |
218 | 218 |
std::vector<IntPair> _arc_vec; |
219 | 219 |
std::vector<Cost> _cost_vec; |
220 | 220 |
IntVector _id_vec; |
221 | 221 |
CostArcMap _cost_map; |
222 | 222 |
CostNodeMap _pi_map; |
223 | 223 |
|
224 | 224 |
public: |
225 | 225 |
|
226 | 226 |
/// \brief Constant for infinite upper bounds (capacities). |
227 | 227 |
/// |
228 | 228 |
/// Constant for infinite upper bounds (capacities). |
229 | 229 |
/// It is \c std::numeric_limits<Value>::infinity() if available, |
230 | 230 |
/// \c std::numeric_limits<Value>::max() otherwise. |
231 | 231 |
const Value INF; |
232 | 232 |
|
233 | 233 |
public: |
234 | 234 |
|
235 | 235 |
/// \brief Constructor. |
236 | 236 |
/// |
237 | 237 |
/// The constructor of the class. |
238 | 238 |
/// |
239 | 239 |
/// \param graph The digraph the algorithm runs on. |
240 | 240 |
CycleCanceling(const GR& graph) : |
241 | 241 |
_graph(graph), _node_id(graph), _arc_idf(graph), _arc_idb(graph), |
242 | 242 |
_cost_map(_cost_vec), _pi_map(_pi), |
243 | 243 |
INF(std::numeric_limits<Value>::has_infinity ? |
244 | 244 |
std::numeric_limits<Value>::infinity() : |
245 | 245 |
std::numeric_limits<Value>::max()) |
246 | 246 |
{ |
247 | 247 |
// Check the number types |
248 | 248 |
LEMON_ASSERT(std::numeric_limits<Value>::is_signed, |
249 | 249 |
"The flow type of CycleCanceling must be signed"); |
250 | 250 |
LEMON_ASSERT(std::numeric_limits<Cost>::is_signed, |
251 | 251 |
"The cost type of CycleCanceling must be signed"); |
252 | 252 |
|
253 |
// Resize vectors |
|
254 |
_node_num = countNodes(_graph); |
|
255 |
_arc_num = countArcs(_graph); |
|
256 |
_res_node_num = _node_num + 1; |
|
257 |
_res_arc_num = 2 * (_arc_num + _node_num); |
|
258 |
_root = _node_num; |
|
259 |
|
|
260 |
_first_out.resize(_res_node_num + 1); |
|
261 |
_forward.resize(_res_arc_num); |
|
262 |
_source.resize(_res_arc_num); |
|
263 |
_target.resize(_res_arc_num); |
|
264 |
_reverse.resize(_res_arc_num); |
|
265 |
|
|
266 |
_lower.resize(_res_arc_num); |
|
267 |
_upper.resize(_res_arc_num); |
|
268 |
_cost.resize(_res_arc_num); |
|
269 |
_supply.resize(_res_node_num); |
|
270 |
|
|
271 |
_res_cap.resize(_res_arc_num); |
|
272 |
_pi.resize(_res_node_num); |
|
273 |
|
|
274 |
_arc_vec.reserve(_res_arc_num); |
|
275 |
_cost_vec.reserve(_res_arc_num); |
|
276 |
_id_vec.reserve(_res_arc_num); |
|
277 |
|
|
278 |
// Copy the graph |
|
279 |
int i = 0, j = 0, k = 2 * _arc_num + _node_num; |
|
280 |
for (NodeIt n(_graph); n != INVALID; ++n, ++i) { |
|
281 |
_node_id[n] = i; |
|
282 |
} |
|
283 |
i = 0; |
|
284 |
for (NodeIt n(_graph); n != INVALID; ++n, ++i) { |
|
285 |
_first_out[i] = j; |
|
286 |
for (OutArcIt a(_graph, n); a != INVALID; ++a, ++j) { |
|
287 |
_arc_idf[a] = j; |
|
288 |
_forward[j] = true; |
|
289 |
_source[j] = i; |
|
290 |
_target[j] = _node_id[_graph.runningNode(a)]; |
|
291 |
} |
|
292 |
for (InArcIt a(_graph, n); a != INVALID; ++a, ++j) { |
|
293 |
_arc_idb[a] = j; |
|
294 |
_forward[j] = false; |
|
295 |
_source[j] = i; |
|
296 |
_target[j] = _node_id[_graph.runningNode(a)]; |
|
297 |
} |
|
298 |
_forward[j] = false; |
|
299 |
_source[j] = i; |
|
300 |
_target[j] = _root; |
|
301 |
_reverse[j] = k; |
|
302 |
_forward[k] = true; |
|
303 |
_source[k] = _root; |
|
304 |
_target[k] = i; |
|
305 |
_reverse[k] = j; |
|
306 |
++j; ++k; |
|
307 |
} |
|
308 |
_first_out[i] = j; |
|
309 |
_first_out[_res_node_num] = k; |
|
310 |
for (ArcIt a(_graph); a != INVALID; ++a) { |
|
311 |
int fi = _arc_idf[a]; |
|
312 |
int bi = _arc_idb[a]; |
|
313 |
_reverse[fi] = bi; |
|
314 |
_reverse[bi] = fi; |
|
315 |
} |
|
316 |
|
|
317 |
// Reset |
|
253 |
// Reset data structures |
|
318 | 254 |
reset(); |
319 | 255 |
} |
320 | 256 |
|
321 | 257 |
/// \name Parameters |
322 | 258 |
/// The parameters of the algorithm can be specified using these |
323 | 259 |
/// functions. |
324 | 260 |
|
325 | 261 |
/// @{ |
326 | 262 |
|
327 | 263 |
/// \brief Set the lower bounds on the arcs. |
328 | 264 |
/// |
329 | 265 |
/// This function sets the lower bounds on the arcs. |
330 | 266 |
/// If it is not used before calling \ref run(), the lower bounds |
331 | 267 |
/// will be set to zero on all arcs. |
332 | 268 |
/// |
333 | 269 |
/// \param map An arc map storing the lower bounds. |
334 | 270 |
/// Its \c Value type must be convertible to the \c Value type |
335 | 271 |
/// of the algorithm. |
336 | 272 |
/// |
337 | 273 |
/// \return <tt>(*this)</tt> |
338 | 274 |
template <typename LowerMap> |
339 | 275 |
CycleCanceling& lowerMap(const LowerMap& map) { |
340 | 276 |
_have_lower = true; |
341 | 277 |
for (ArcIt a(_graph); a != INVALID; ++a) { |
342 | 278 |
_lower[_arc_idf[a]] = map[a]; |
343 | 279 |
_lower[_arc_idb[a]] = map[a]; |
344 | 280 |
} |
345 | 281 |
return *this; |
346 | 282 |
} |
347 | 283 |
|
348 | 284 |
/// \brief Set the upper bounds (capacities) on the arcs. |
349 | 285 |
/// |
350 | 286 |
/// This function sets the upper bounds (capacities) on the arcs. |
351 | 287 |
/// If it is not used before calling \ref run(), the upper bounds |
352 | 288 |
/// will be set to \ref INF on all arcs (i.e. the flow value will be |
353 | 289 |
/// unbounded from above). |
354 | 290 |
/// |
355 | 291 |
/// \param map An arc map storing the upper bounds. |
356 | 292 |
/// Its \c Value type must be convertible to the \c Value type |
357 | 293 |
/// of the algorithm. |
358 | 294 |
/// |
359 | 295 |
/// \return <tt>(*this)</tt> |
360 | 296 |
template<typename UpperMap> |
361 | 297 |
CycleCanceling& upperMap(const UpperMap& map) { |
362 | 298 |
for (ArcIt a(_graph); a != INVALID; ++a) { |
363 | 299 |
_upper[_arc_idf[a]] = map[a]; |
364 | 300 |
} |
365 | 301 |
return *this; |
366 | 302 |
} |
367 | 303 |
|
368 | 304 |
/// \brief Set the costs of the arcs. |
369 | 305 |
/// |
370 | 306 |
/// This function sets the costs of the arcs. |
371 | 307 |
/// If it is not used before calling \ref run(), the costs |
372 | 308 |
/// will be set to \c 1 on all arcs. |
373 | 309 |
/// |
374 | 310 |
/// \param map An arc map storing the costs. |
375 | 311 |
/// Its \c Value type must be convertible to the \c Cost type |
376 | 312 |
/// of the algorithm. |
377 | 313 |
/// |
378 | 314 |
/// \return <tt>(*this)</tt> |
379 | 315 |
template<typename CostMap> |
380 | 316 |
CycleCanceling& costMap(const CostMap& map) { |
381 | 317 |
for (ArcIt a(_graph); a != INVALID; ++a) { |
382 | 318 |
_cost[_arc_idf[a]] = map[a]; |
383 | 319 |
_cost[_arc_idb[a]] = -map[a]; |
384 | 320 |
} |
385 | 321 |
return *this; |
386 | 322 |
} |
387 | 323 |
|
388 | 324 |
/// \brief Set the supply values of the nodes. |
389 | 325 |
/// |
390 | 326 |
/// This function sets the supply values of the nodes. |
391 | 327 |
/// If neither this function nor \ref stSupply() is used before |
392 | 328 |
/// calling \ref run(), the supply of each node will be set to zero. |
393 | 329 |
/// |
394 | 330 |
/// \param map A node map storing the supply values. |
395 | 331 |
/// Its \c Value type must be convertible to the \c Value type |
396 | 332 |
/// of the algorithm. |
397 | 333 |
/// |
398 | 334 |
/// \return <tt>(*this)</tt> |
399 | 335 |
template<typename SupplyMap> |
400 | 336 |
CycleCanceling& supplyMap(const SupplyMap& map) { |
401 | 337 |
for (NodeIt n(_graph); n != INVALID; ++n) { |
402 | 338 |
_supply[_node_id[n]] = map[n]; |
403 | 339 |
} |
404 | 340 |
return *this; |
405 | 341 |
} |
406 | 342 |
|
407 | 343 |
/// \brief Set single source and target nodes and a supply value. |
408 | 344 |
/// |
409 | 345 |
/// This function sets a single source node and a single target node |
410 | 346 |
/// and the required flow value. |
411 | 347 |
/// If neither this function nor \ref supplyMap() is used before |
412 | 348 |
/// calling \ref run(), the supply of each node will be set to zero. |
413 | 349 |
/// |
414 | 350 |
/// Using this function has the same effect as using \ref supplyMap() |
415 | 351 |
/// with such a map in which \c k is assigned to \c s, \c -k is |
416 | 352 |
/// assigned to \c t and all other nodes have zero supply value. |
417 | 353 |
/// |
418 | 354 |
/// \param s The source node. |
419 | 355 |
/// \param t The target node. |
420 | 356 |
/// \param k The required amount of flow from node \c s to node \c t |
421 | 357 |
/// (i.e. the supply of \c s and the demand of \c t). |
422 | 358 |
/// |
423 | 359 |
/// \return <tt>(*this)</tt> |
424 | 360 |
CycleCanceling& stSupply(const Node& s, const Node& t, Value k) { |
425 | 361 |
for (int i = 0; i != _res_node_num; ++i) { |
426 | 362 |
_supply[i] = 0; |
427 | 363 |
} |
428 | 364 |
_supply[_node_id[s]] = k; |
429 | 365 |
_supply[_node_id[t]] = -k; |
430 | 366 |
return *this; |
431 | 367 |
} |
432 | 368 |
|
433 | 369 |
/// @} |
434 | 370 |
|
435 | 371 |
/// \name Execution control |
436 | 372 |
/// The algorithm can be executed using \ref run(). |
437 | 373 |
|
438 | 374 |
/// @{ |
439 | 375 |
|
440 | 376 |
/// \brief Run the algorithm. |
441 | 377 |
/// |
442 | 378 |
/// This function runs the algorithm. |
443 | 379 |
/// The paramters can be specified using functions \ref lowerMap(), |
444 | 380 |
/// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply(). |
445 | 381 |
/// For example, |
446 | 382 |
/// \code |
447 | 383 |
/// CycleCanceling<ListDigraph> cc(graph); |
448 | 384 |
/// cc.lowerMap(lower).upperMap(upper).costMap(cost) |
449 | 385 |
/// .supplyMap(sup).run(); |
450 | 386 |
/// \endcode |
451 | 387 |
/// |
452 |
/// This function can be called more than once. All the parameters |
|
453 |
/// that have been given are kept for the next call, unless |
|
454 |
/// \ref reset() is called, thus only the modified parameters |
|
455 |
/// have to be set again. See \ref reset() for examples. |
|
456 |
/// However, the underlying digraph must not be modified after this |
|
457 |
/// class have been constructed, since it copies and extends the graph. |
|
388 |
/// This function can be called more than once. All the given parameters |
|
389 |
/// are kept for the next call, unless \ref resetParams() or \ref reset() |
|
390 |
/// is used, thus only the modified parameters have to be set again. |
|
391 |
/// If the underlying digraph was also modified after the construction |
|
392 |
/// of the class (or the last \ref reset() call), then the \ref reset() |
|
393 |
/// function must be called. |
|
458 | 394 |
/// |
459 | 395 |
/// \param method The cycle-canceling method that will be used. |
460 | 396 |
/// For more information, see \ref Method. |
461 | 397 |
/// |
462 | 398 |
/// \return \c INFEASIBLE if no feasible flow exists, |
463 | 399 |
/// \n \c OPTIMAL if the problem has optimal solution |
464 | 400 |
/// (i.e. it is feasible and bounded), and the algorithm has found |
465 | 401 |
/// optimal flow and node potentials (primal and dual solutions), |
466 | 402 |
/// \n \c UNBOUNDED if the digraph contains an arc of negative cost |
467 | 403 |
/// and infinite upper bound. It means that the objective function |
468 | 404 |
/// is unbounded on that arc, however, note that it could actually be |
469 | 405 |
/// bounded over the feasible flows, but this algroithm cannot handle |
470 | 406 |
/// these cases. |
471 | 407 |
/// |
472 | 408 |
/// \see ProblemType, Method |
409 |
/// \see resetParams(), reset() |
|
473 | 410 |
ProblemType run(Method method = CANCEL_AND_TIGHTEN) { |
474 | 411 |
ProblemType pt = init(); |
475 | 412 |
if (pt != OPTIMAL) return pt; |
476 | 413 |
start(method); |
477 | 414 |
return OPTIMAL; |
478 | 415 |
} |
479 | 416 |
|
480 | 417 |
/// \brief Reset all the parameters that have been given before. |
481 | 418 |
/// |
482 | 419 |
/// This function resets all the paramaters that have been given |
483 | 420 |
/// before using functions \ref lowerMap(), \ref upperMap(), |
484 | 421 |
/// \ref costMap(), \ref supplyMap(), \ref stSupply(). |
485 | 422 |
/// |
486 |
/// It is useful for multiple run() calls. If this function is not |
|
487 |
/// used, all the parameters given before are kept for the next |
|
488 |
/// \ref run() call. |
|
489 |
/// However, the underlying digraph must not be modified after this |
|
490 |
/// |
|
423 |
/// It is useful for multiple \ref run() calls. Basically, all the given |
|
424 |
/// parameters are kept for the next \ref run() call, unless |
|
425 |
/// \ref resetParams() or \ref reset() is used. |
|
426 |
/// If the underlying digraph was also modified after the construction |
|
427 |
/// of the class or the last \ref reset() call, then the \ref reset() |
|
428 |
/// function must be used, otherwise \ref resetParams() is sufficient. |
|
491 | 429 |
/// |
492 | 430 |
/// For example, |
493 | 431 |
/// \code |
494 | 432 |
/// CycleCanceling<ListDigraph> cs(graph); |
495 | 433 |
/// |
496 | 434 |
/// // First run |
497 | 435 |
/// cc.lowerMap(lower).upperMap(upper).costMap(cost) |
498 | 436 |
/// .supplyMap(sup).run(); |
499 | 437 |
/// |
500 |
/// // Run again with modified cost map ( |
|
438 |
/// // Run again with modified cost map (resetParams() is not called, |
|
501 | 439 |
/// // so only the cost map have to be set again) |
502 | 440 |
/// cost[e] += 100; |
503 | 441 |
/// cc.costMap(cost).run(); |
504 | 442 |
/// |
505 |
/// // Run again from scratch using |
|
443 |
/// // Run again from scratch using resetParams() |
|
506 | 444 |
/// // (the lower bounds will be set to zero on all arcs) |
507 |
/// cc. |
|
445 |
/// cc.resetParams(); |
|
508 | 446 |
/// cc.upperMap(capacity).costMap(cost) |
509 | 447 |
/// .supplyMap(sup).run(); |
510 | 448 |
/// \endcode |
511 | 449 |
/// |
512 | 450 |
/// \return <tt>(*this)</tt> |
513 |
|
|
451 |
/// |
|
452 |
/// \see reset(), run() |
|
453 |
CycleCanceling& resetParams() { |
|
514 | 454 |
for (int i = 0; i != _res_node_num; ++i) { |
515 | 455 |
_supply[i] = 0; |
516 | 456 |
} |
517 | 457 |
int limit = _first_out[_root]; |
518 | 458 |
for (int j = 0; j != limit; ++j) { |
519 | 459 |
_lower[j] = 0; |
520 | 460 |
_upper[j] = INF; |
521 | 461 |
_cost[j] = _forward[j] ? 1 : -1; |
522 | 462 |
} |
523 | 463 |
for (int j = limit; j != _res_arc_num; ++j) { |
524 | 464 |
_lower[j] = 0; |
525 | 465 |
_upper[j] = INF; |
526 | 466 |
_cost[j] = 0; |
527 | 467 |
_cost[_reverse[j]] = 0; |
528 | 468 |
} |
529 | 469 |
_have_lower = false; |
530 | 470 |
return *this; |
531 | 471 |
} |
532 | 472 |
|
473 |
/// \brief Reset the internal data structures and all the parameters |
|
474 |
/// that have been given before. |
|
475 |
/// |
|
476 |
/// This function resets the internal data structures and all the |
|
477 |
/// paramaters that have been given before using functions \ref lowerMap(), |
|
478 |
/// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply(). |
|
479 |
/// |
|
480 |
/// It is useful for multiple \ref run() calls. Basically, all the given |
|
481 |
/// parameters are kept for the next \ref run() call, unless |
|
482 |
/// \ref resetParams() or \ref reset() is used. |
|
483 |
/// If the underlying digraph was also modified after the construction |
|
484 |
/// of the class or the last \ref reset() call, then the \ref reset() |
|
485 |
/// function must be used, otherwise \ref resetParams() is sufficient. |
|
486 |
/// |
|
487 |
/// See \ref resetParams() for examples. |
|
488 |
/// |
|
489 |
/// \return <tt>(*this)</tt> |
|
490 |
/// |
|
491 |
/// \see resetParams(), run() |
|
492 |
CycleCanceling& reset() { |
|
493 |
// Resize vectors |
|
494 |
_node_num = countNodes(_graph); |
|
495 |
_arc_num = countArcs(_graph); |
|
496 |
_res_node_num = _node_num + 1; |
|
497 |
_res_arc_num = 2 * (_arc_num + _node_num); |
|
498 |
_root = _node_num; |
|
499 |
|
|
500 |
_first_out.resize(_res_node_num + 1); |
|
501 |
_forward.resize(_res_arc_num); |
|
502 |
_source.resize(_res_arc_num); |
|
503 |
_target.resize(_res_arc_num); |
|
504 |
_reverse.resize(_res_arc_num); |
|
505 |
|
|
506 |
_lower.resize(_res_arc_num); |
|
507 |
_upper.resize(_res_arc_num); |
|
508 |
_cost.resize(_res_arc_num); |
|
509 |
_supply.resize(_res_node_num); |
|
510 |
|
|
511 |
_res_cap.resize(_res_arc_num); |
|
512 |
_pi.resize(_res_node_num); |
|
513 |
|
|
514 |
_arc_vec.reserve(_res_arc_num); |
|
515 |
_cost_vec.reserve(_res_arc_num); |
|
516 |
_id_vec.reserve(_res_arc_num); |
|
517 |
|
|
518 |
// Copy the graph |
|
519 |
int i = 0, j = 0, k = 2 * _arc_num + _node_num; |
|
520 |
for (NodeIt n(_graph); n != INVALID; ++n, ++i) { |
|
521 |
_node_id[n] = i; |
|
522 |
} |
|
523 |
i = 0; |
|
524 |
for (NodeIt n(_graph); n != INVALID; ++n, ++i) { |
|
525 |
_first_out[i] = j; |
|
526 |
for (OutArcIt a(_graph, n); a != INVALID; ++a, ++j) { |
|
527 |
_arc_idf[a] = j; |
|
528 |
_forward[j] = true; |
|
529 |
_source[j] = i; |
|
530 |
_target[j] = _node_id[_graph.runningNode(a)]; |
|
531 |
} |
|
532 |
for (InArcIt a(_graph, n); a != INVALID; ++a, ++j) { |
|
533 |
_arc_idb[a] = j; |
|
534 |
_forward[j] = false; |
|
535 |
_source[j] = i; |
|
536 |
_target[j] = _node_id[_graph.runningNode(a)]; |
|
537 |
} |
|
538 |
_forward[j] = false; |
|
539 |
_source[j] = i; |
|
540 |
_target[j] = _root; |
|
541 |
_reverse[j] = k; |
|
542 |
_forward[k] = true; |
|
543 |
_source[k] = _root; |
|
544 |
_target[k] = i; |
|
545 |
_reverse[k] = j; |
|
546 |
++j; ++k; |
|
547 |
} |
|
548 |
_first_out[i] = j; |
|
549 |
_first_out[_res_node_num] = k; |
|
550 |
for (ArcIt a(_graph); a != INVALID; ++a) { |
|
551 |
int fi = _arc_idf[a]; |
|
552 |
int bi = _arc_idb[a]; |
|
553 |
_reverse[fi] = bi; |
|
554 |
_reverse[bi] = fi; |
|
555 |
} |
|
556 |
|
|
557 |
// Reset parameters |
|
558 |
resetParams(); |
|
559 |
return *this; |
|
560 |
} |
|
561 |
|
|
533 | 562 |
/// @} |
534 | 563 |
|
535 | 564 |
/// \name Query Functions |
536 | 565 |
/// The results of the algorithm can be obtained using these |
537 | 566 |
/// functions.\n |
538 | 567 |
/// The \ref run() function must be called before using them. |
539 | 568 |
|
540 | 569 |
/// @{ |
541 | 570 |
|
542 | 571 |
/// \brief Return the total cost of the found flow. |
543 | 572 |
/// |
544 | 573 |
/// This function returns the total cost of the found flow. |
545 | 574 |
/// Its complexity is O(e). |
546 | 575 |
/// |
547 | 576 |
/// \note The return type of the function can be specified as a |
548 | 577 |
/// template parameter. For example, |
549 | 578 |
/// \code |
550 | 579 |
/// cc.totalCost<double>(); |
551 | 580 |
/// \endcode |
552 | 581 |
/// It is useful if the total cost cannot be stored in the \c Cost |
553 | 582 |
/// type of the algorithm, which is the default return type of the |
554 | 583 |
/// function. |
555 | 584 |
/// |
556 | 585 |
/// \pre \ref run() must be called before using this function. |
557 | 586 |
template <typename Number> |
558 | 587 |
Number totalCost() const { |
559 | 588 |
Number c = 0; |
560 | 589 |
for (ArcIt a(_graph); a != INVALID; ++a) { |
561 | 590 |
int i = _arc_idb[a]; |
562 | 591 |
c += static_cast<Number>(_res_cap[i]) * |
563 | 592 |
(-static_cast<Number>(_cost[i])); |
564 | 593 |
} |
565 | 594 |
return c; |
566 | 595 |
} |
567 | 596 |
|
568 | 597 |
#ifndef DOXYGEN |
569 | 598 |
Cost totalCost() const { |
570 | 599 |
return totalCost<Cost>(); |
571 | 600 |
} |
572 | 601 |
#endif |
573 | 602 |
|
574 | 603 |
/// \brief Return the flow on the given arc. |
575 | 604 |
/// |
576 | 605 |
/// This function returns the flow on the given arc. |
577 | 606 |
/// |
578 | 607 |
/// \pre \ref run() must be called before using this function. |
579 | 608 |
Value flow(const Arc& a) const { |
580 | 609 |
return _res_cap[_arc_idb[a]]; |
581 | 610 |
} |
582 | 611 |
|
583 | 612 |
/// \brief Return the flow map (the primal solution). |
584 | 613 |
/// |
585 | 614 |
/// This function copies the flow value on each arc into the given |
586 | 615 |
/// map. The \c Value type of the algorithm must be convertible to |
587 | 616 |
/// the \c Value type of the map. |
588 | 617 |
/// |
589 | 618 |
/// \pre \ref run() must be called before using this function. |
590 | 619 |
template <typename FlowMap> |
591 | 620 |
void flowMap(FlowMap &map) const { |
592 | 621 |
for (ArcIt a(_graph); a != INVALID; ++a) { |
593 | 622 |
map.set(a, _res_cap[_arc_idb[a]]); |
594 | 623 |
} |
595 | 624 |
} |
596 | 625 |
|
597 | 626 |
/// \brief Return the potential (dual value) of the given node. |
598 | 627 |
/// |
599 | 628 |
/// This function returns the potential (dual value) of the |
600 | 629 |
/// given node. |
601 | 630 |
/// |
602 | 631 |
/// \pre \ref run() must be called before using this function. |
603 | 632 |
Cost potential(const Node& n) const { |
604 | 633 |
return static_cast<Cost>(_pi[_node_id[n]]); |
605 | 634 |
} |
606 | 635 |
|
607 | 636 |
/// \brief Return the potential map (the dual solution). |
608 | 637 |
/// |
609 | 638 |
/// This function copies the potential (dual value) of each node |
610 | 639 |
/// into the given map. |
611 | 640 |
/// The \c Cost type of the algorithm must be convertible to the |
612 | 641 |
/// \c Value type of the map. |
613 | 642 |
/// |
614 | 643 |
/// \pre \ref run() must be called before using this function. |
615 | 644 |
template <typename PotentialMap> |
616 | 645 |
void potentialMap(PotentialMap &map) const { |
617 | 646 |
for (NodeIt n(_graph); n != INVALID; ++n) { |
618 | 647 |
map.set(n, static_cast<Cost>(_pi[_node_id[n]])); |
619 | 648 |
} |
620 | 649 |
} |
621 | 650 |
|
622 | 651 |
/// @} |
623 | 652 |
|
624 | 653 |
private: |
625 | 654 |
|
626 | 655 |
// Initialize the algorithm |
627 | 656 |
ProblemType init() { |
628 | 657 |
if (_res_node_num <= 1) return INFEASIBLE; |
... | ... |
@@ -101,192 +101,193 @@ |
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 | 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 |
bool _arc_mixing; |
|
197 | 198 |
|
198 | 199 |
// Node and arc data |
199 | 200 |
ValueVector _lower; |
200 | 201 |
ValueVector _upper; |
201 | 202 |
ValueVector _cap; |
202 | 203 |
CostVector _cost; |
203 | 204 |
ValueVector _supply; |
204 | 205 |
ValueVector _flow; |
205 | 206 |
CostVector _pi; |
206 | 207 |
|
207 | 208 |
// Data for storing the spanning tree structure |
208 | 209 |
IntVector _parent; |
209 | 210 |
IntVector _pred; |
210 | 211 |
IntVector _thread; |
211 | 212 |
IntVector _rev_thread; |
212 | 213 |
IntVector _succ_num; |
213 | 214 |
IntVector _last_succ; |
214 | 215 |
IntVector _dirty_revs; |
215 | 216 |
CharVector _forward; |
216 | 217 |
CharVector _state; |
217 | 218 |
int _root; |
218 | 219 |
|
219 | 220 |
// Temporary data used in the current pivot iteration |
220 | 221 |
int in_arc, join, u_in, v_in, u_out, v_out; |
221 | 222 |
int first, second, right, last; |
222 | 223 |
int stem, par_stem, new_stem; |
223 | 224 |
Value delta; |
224 | 225 |
|
225 | 226 |
const Value MAX; |
226 | 227 |
|
227 | 228 |
public: |
228 | 229 |
|
229 | 230 |
/// \brief Constant for infinite upper bounds (capacities). |
230 | 231 |
/// |
231 | 232 |
/// Constant for infinite upper bounds (capacities). |
232 | 233 |
/// It is \c std::numeric_limits<Value>::infinity() if available, |
233 | 234 |
/// \c std::numeric_limits<Value>::max() otherwise. |
234 | 235 |
const Value INF; |
235 | 236 |
|
236 | 237 |
private: |
237 | 238 |
|
238 | 239 |
// Implementation of the First Eligible pivot rule |
239 | 240 |
class FirstEligiblePivotRule |
240 | 241 |
{ |
241 | 242 |
private: |
242 | 243 |
|
243 | 244 |
// References to the NetworkSimplex class |
244 | 245 |
const IntVector &_source; |
245 | 246 |
const IntVector &_target; |
246 | 247 |
const CostVector &_cost; |
247 | 248 |
const CharVector &_state; |
248 | 249 |
const CostVector &_pi; |
249 | 250 |
int &_in_arc; |
250 | 251 |
int _search_arc_num; |
251 | 252 |
|
252 | 253 |
// Pivot rule data |
253 | 254 |
int _next_arc; |
254 | 255 |
|
255 | 256 |
public: |
256 | 257 |
|
257 | 258 |
// Constructor |
258 | 259 |
FirstEligiblePivotRule(NetworkSimplex &ns) : |
259 | 260 |
_source(ns._source), _target(ns._target), |
260 | 261 |
_cost(ns._cost), _state(ns._state), _pi(ns._pi), |
261 | 262 |
_in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num), |
262 | 263 |
_next_arc(0) |
263 | 264 |
{} |
264 | 265 |
|
265 | 266 |
// Find next entering arc |
266 | 267 |
bool findEnteringArc() { |
267 | 268 |
Cost c; |
268 | 269 |
for (int e = _next_arc; e < _search_arc_num; ++e) { |
269 | 270 |
c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
270 | 271 |
if (c < 0) { |
271 | 272 |
_in_arc = e; |
272 | 273 |
_next_arc = e + 1; |
273 | 274 |
return true; |
274 | 275 |
} |
275 | 276 |
} |
276 | 277 |
for (int e = 0; e < _next_arc; ++e) { |
277 | 278 |
c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
278 | 279 |
if (c < 0) { |
279 | 280 |
_in_arc = e; |
280 | 281 |
_next_arc = e + 1; |
281 | 282 |
return true; |
282 | 283 |
} |
283 | 284 |
} |
284 | 285 |
return false; |
285 | 286 |
} |
286 | 287 |
|
287 | 288 |
}; //class FirstEligiblePivotRule |
288 | 289 |
|
289 | 290 |
|
290 | 291 |
// Implementation of the Best Eligible pivot rule |
291 | 292 |
class BestEligiblePivotRule |
292 | 293 |
{ |
... | ... |
@@ -540,472 +541,503 @@ |
540 | 541 |
|
541 | 542 |
// Constructor |
542 | 543 |
AlteringListPivotRule(NetworkSimplex &ns) : |
543 | 544 |
_source(ns._source), _target(ns._target), |
544 | 545 |
_cost(ns._cost), _state(ns._state), _pi(ns._pi), |
545 | 546 |
_in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num), |
546 | 547 |
_next_arc(0), _cand_cost(ns._search_arc_num), _sort_func(_cand_cost) |
547 | 548 |
{ |
548 | 549 |
// The main parameters of the pivot rule |
549 | 550 |
const double BLOCK_SIZE_FACTOR = 1.0; |
550 | 551 |
const int MIN_BLOCK_SIZE = 10; |
551 | 552 |
const double HEAD_LENGTH_FACTOR = 0.1; |
552 | 553 |
const int MIN_HEAD_LENGTH = 3; |
553 | 554 |
|
554 | 555 |
_block_size = std::max( int(BLOCK_SIZE_FACTOR * |
555 | 556 |
std::sqrt(double(_search_arc_num))), |
556 | 557 |
MIN_BLOCK_SIZE ); |
557 | 558 |
_head_length = std::max( int(HEAD_LENGTH_FACTOR * _block_size), |
558 | 559 |
MIN_HEAD_LENGTH ); |
559 | 560 |
_candidates.resize(_head_length + _block_size); |
560 | 561 |
_curr_length = 0; |
561 | 562 |
} |
562 | 563 |
|
563 | 564 |
// Find next entering arc |
564 | 565 |
bool findEnteringArc() { |
565 | 566 |
// Check the current candidate list |
566 | 567 |
int e; |
567 | 568 |
for (int i = 0; i < _curr_length; ++i) { |
568 | 569 |
e = _candidates[i]; |
569 | 570 |
_cand_cost[e] = _state[e] * |
570 | 571 |
(_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
571 | 572 |
if (_cand_cost[e] >= 0) { |
572 | 573 |
_candidates[i--] = _candidates[--_curr_length]; |
573 | 574 |
} |
574 | 575 |
} |
575 | 576 |
|
576 | 577 |
// Extend the list |
577 | 578 |
int cnt = _block_size; |
578 | 579 |
int limit = _head_length; |
579 | 580 |
|
580 | 581 |
for (e = _next_arc; e < _search_arc_num; ++e) { |
581 | 582 |
_cand_cost[e] = _state[e] * |
582 | 583 |
(_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
583 | 584 |
if (_cand_cost[e] < 0) { |
584 | 585 |
_candidates[_curr_length++] = e; |
585 | 586 |
} |
586 | 587 |
if (--cnt == 0) { |
587 | 588 |
if (_curr_length > limit) goto search_end; |
588 | 589 |
limit = 0; |
589 | 590 |
cnt = _block_size; |
590 | 591 |
} |
591 | 592 |
} |
592 | 593 |
for (e = 0; e < _next_arc; ++e) { |
593 | 594 |
_cand_cost[e] = _state[e] * |
594 | 595 |
(_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
595 | 596 |
if (_cand_cost[e] < 0) { |
596 | 597 |
_candidates[_curr_length++] = e; |
597 | 598 |
} |
598 | 599 |
if (--cnt == 0) { |
599 | 600 |
if (_curr_length > limit) goto search_end; |
600 | 601 |
limit = 0; |
601 | 602 |
cnt = _block_size; |
602 | 603 |
} |
603 | 604 |
} |
604 | 605 |
if (_curr_length == 0) return false; |
605 | 606 |
|
606 | 607 |
search_end: |
607 | 608 |
|
608 | 609 |
// Make heap of the candidate list (approximating a partial sort) |
609 | 610 |
make_heap( _candidates.begin(), _candidates.begin() + _curr_length, |
610 | 611 |
_sort_func ); |
611 | 612 |
|
612 | 613 |
// Pop the first element of the heap |
613 | 614 |
_in_arc = _candidates[0]; |
614 | 615 |
_next_arc = e; |
615 | 616 |
pop_heap( _candidates.begin(), _candidates.begin() + _curr_length, |
616 | 617 |
_sort_func ); |
617 | 618 |
_curr_length = std::min(_head_length, _curr_length - 1); |
618 | 619 |
return true; |
619 | 620 |
} |
620 | 621 |
|
621 | 622 |
}; //class AlteringListPivotRule |
622 | 623 |
|
623 | 624 |
public: |
624 | 625 |
|
625 | 626 |
/// \brief Constructor. |
626 | 627 |
/// |
627 | 628 |
/// The constructor of the class. |
628 | 629 |
/// |
629 | 630 |
/// \param graph The digraph the algorithm runs on. |
630 | 631 |
/// \param arc_mixing Indicate if the arcs have to be stored in a |
631 | 632 |
/// mixed order in the internal data structure. |
632 | 633 |
/// In special cases, it could lead to better overall performance, |
633 | 634 |
/// but it is usually slower. Therefore it is disabled by default. |
634 | 635 |
NetworkSimplex(const GR& graph, bool arc_mixing = false) : |
635 | 636 |
_graph(graph), _node_id(graph), _arc_id(graph), |
637 |
_arc_mixing(arc_mixing), |
|
636 | 638 |
MAX(std::numeric_limits<Value>::max()), |
637 | 639 |
INF(std::numeric_limits<Value>::has_infinity ? |
638 | 640 |
std::numeric_limits<Value>::infinity() : MAX) |
639 | 641 |
{ |
640 | 642 |
// Check the number types |
641 | 643 |
LEMON_ASSERT(std::numeric_limits<Value>::is_signed, |
642 | 644 |
"The flow type of NetworkSimplex must be signed"); |
643 | 645 |
LEMON_ASSERT(std::numeric_limits<Cost>::is_signed, |
644 | 646 |
"The cost type of NetworkSimplex must be signed"); |
645 | 647 |
|
646 |
// Resize vectors |
|
647 |
_node_num = countNodes(_graph); |
|
648 |
_arc_num = countArcs(_graph); |
|
649 |
int all_node_num = _node_num + 1; |
|
650 |
int max_arc_num = _arc_num + 2 * _node_num; |
|
651 |
|
|
652 |
_source.resize(max_arc_num); |
|
653 |
_target.resize(max_arc_num); |
|
654 |
|
|
655 |
_lower.resize(_arc_num); |
|
656 |
_upper.resize(_arc_num); |
|
657 |
_cap.resize(max_arc_num); |
|
658 |
_cost.resize(max_arc_num); |
|
659 |
_supply.resize(all_node_num); |
|
660 |
_flow.resize(max_arc_num); |
|
661 |
_pi.resize(all_node_num); |
|
662 |
|
|
663 |
_parent.resize(all_node_num); |
|
664 |
_pred.resize(all_node_num); |
|
665 |
_forward.resize(all_node_num); |
|
666 |
_thread.resize(all_node_num); |
|
667 |
_rev_thread.resize(all_node_num); |
|
668 |
_succ_num.resize(all_node_num); |
|
669 |
_last_succ.resize(all_node_num); |
|
670 |
_state.resize(max_arc_num); |
|
671 |
|
|
672 |
// Copy the graph |
|
673 |
int i = 0; |
|
674 |
for (NodeIt n(_graph); n != INVALID; ++n, ++i) { |
|
675 |
_node_id[n] = i; |
|
676 |
} |
|
677 |
if (arc_mixing) { |
|
678 |
// Store the arcs in a mixed order |
|
679 |
int k = std::max(int(std::sqrt(double(_arc_num))), 10); |
|
680 |
int i = 0, j = 0; |
|
681 |
for (ArcIt a(_graph); a != INVALID; ++a) { |
|
682 |
_arc_id[a] = i; |
|
683 |
_source[i] = _node_id[_graph.source(a)]; |
|
684 |
_target[i] = _node_id[_graph.target(a)]; |
|
685 |
if ((i += k) >= _arc_num) i = ++j; |
|
686 |
} |
|
687 |
} else { |
|
688 |
// Store the arcs in the original order |
|
689 |
int i = 0; |
|
690 |
for (ArcIt a(_graph); a != INVALID; ++a, ++i) { |
|
691 |
_arc_id[a] = i; |
|
692 |
_source[i] = _node_id[_graph.source(a)]; |
|
693 |
_target[i] = _node_id[_graph.target(a)]; |
|
694 |
} |
|
695 |
} |
|
696 |
|
|
697 |
// Reset parameters |
|
648 |
// Reset data structures |
|
698 | 649 |
reset(); |
699 | 650 |
} |
700 | 651 |
|
701 | 652 |
/// \name Parameters |
702 | 653 |
/// The parameters of the algorithm can be specified using these |
703 | 654 |
/// functions. |
704 | 655 |
|
705 | 656 |
/// @{ |
706 | 657 |
|
707 | 658 |
/// \brief Set the lower bounds on the arcs. |
708 | 659 |
/// |
709 | 660 |
/// This function sets the lower bounds on the arcs. |
710 | 661 |
/// If it is not used before calling \ref run(), the lower bounds |
711 | 662 |
/// will be set to zero on all arcs. |
712 | 663 |
/// |
713 | 664 |
/// \param map An arc map storing the lower bounds. |
714 | 665 |
/// Its \c Value type must be convertible to the \c Value type |
715 | 666 |
/// of the algorithm. |
716 | 667 |
/// |
717 | 668 |
/// \return <tt>(*this)</tt> |
718 | 669 |
template <typename LowerMap> |
719 | 670 |
NetworkSimplex& lowerMap(const LowerMap& map) { |
720 | 671 |
_have_lower = true; |
721 | 672 |
for (ArcIt a(_graph); a != INVALID; ++a) { |
722 | 673 |
_lower[_arc_id[a]] = map[a]; |
723 | 674 |
} |
724 | 675 |
return *this; |
725 | 676 |
} |
726 | 677 |
|
727 | 678 |
/// \brief Set the upper bounds (capacities) on the arcs. |
728 | 679 |
/// |
729 | 680 |
/// This function sets the upper bounds (capacities) on the arcs. |
730 | 681 |
/// If it is not used before calling \ref run(), the upper bounds |
731 | 682 |
/// will be set to \ref INF on all arcs (i.e. the flow value will be |
732 | 683 |
/// unbounded from above). |
733 | 684 |
/// |
734 | 685 |
/// \param map An arc map storing the upper bounds. |
735 | 686 |
/// Its \c Value type must be convertible to the \c Value type |
736 | 687 |
/// of the algorithm. |
737 | 688 |
/// |
738 | 689 |
/// \return <tt>(*this)</tt> |
739 | 690 |
template<typename UpperMap> |
740 | 691 |
NetworkSimplex& upperMap(const UpperMap& map) { |
741 | 692 |
for (ArcIt a(_graph); a != INVALID; ++a) { |
742 | 693 |
_upper[_arc_id[a]] = map[a]; |
743 | 694 |
} |
744 | 695 |
return *this; |
745 | 696 |
} |
746 | 697 |
|
747 | 698 |
/// \brief Set the costs of the arcs. |
748 | 699 |
/// |
749 | 700 |
/// This function sets the costs of the arcs. |
750 | 701 |
/// If it is not used before calling \ref run(), the costs |
751 | 702 |
/// will be set to \c 1 on all arcs. |
752 | 703 |
/// |
753 | 704 |
/// \param map An arc map storing the costs. |
754 | 705 |
/// Its \c Value type must be convertible to the \c Cost type |
755 | 706 |
/// of the algorithm. |
756 | 707 |
/// |
757 | 708 |
/// \return <tt>(*this)</tt> |
758 | 709 |
template<typename CostMap> |
759 | 710 |
NetworkSimplex& costMap(const CostMap& map) { |
760 | 711 |
for (ArcIt a(_graph); a != INVALID; ++a) { |
761 | 712 |
_cost[_arc_id[a]] = map[a]; |
762 | 713 |
} |
763 | 714 |
return *this; |
764 | 715 |
} |
765 | 716 |
|
766 | 717 |
/// \brief Set the supply values of the nodes. |
767 | 718 |
/// |
768 | 719 |
/// This function sets the supply values of the nodes. |
769 | 720 |
/// If neither this function nor \ref stSupply() is used before |
770 | 721 |
/// calling \ref run(), the supply of each node will be set to zero. |
771 | 722 |
/// |
772 | 723 |
/// \param map A node map storing the supply values. |
773 | 724 |
/// Its \c Value type must be convertible to the \c Value type |
774 | 725 |
/// of the algorithm. |
775 | 726 |
/// |
776 | 727 |
/// \return <tt>(*this)</tt> |
777 | 728 |
template<typename SupplyMap> |
778 | 729 |
NetworkSimplex& supplyMap(const SupplyMap& map) { |
779 | 730 |
for (NodeIt n(_graph); n != INVALID; ++n) { |
780 | 731 |
_supply[_node_id[n]] = map[n]; |
781 | 732 |
} |
782 | 733 |
return *this; |
783 | 734 |
} |
784 | 735 |
|
785 | 736 |
/// \brief Set single source and target nodes and a supply value. |
786 | 737 |
/// |
787 | 738 |
/// This function sets a single source node and a single target node |
788 | 739 |
/// and the required flow value. |
789 | 740 |
/// If neither this function nor \ref supplyMap() is used before |
790 | 741 |
/// calling \ref run(), the supply of each node will be set to zero. |
791 | 742 |
/// |
792 | 743 |
/// Using this function has the same effect as using \ref supplyMap() |
793 | 744 |
/// with such a map in which \c k is assigned to \c s, \c -k is |
794 | 745 |
/// assigned to \c t and all other nodes have zero supply value. |
795 | 746 |
/// |
796 | 747 |
/// \param s The source node. |
797 | 748 |
/// \param t The target node. |
798 | 749 |
/// \param k The required amount of flow from node \c s to node \c t |
799 | 750 |
/// (i.e. the supply of \c s and the demand of \c t). |
800 | 751 |
/// |
801 | 752 |
/// \return <tt>(*this)</tt> |
802 | 753 |
NetworkSimplex& stSupply(const Node& s, const Node& t, Value k) { |
803 | 754 |
for (int i = 0; i != _node_num; ++i) { |
804 | 755 |
_supply[i] = 0; |
805 | 756 |
} |
806 | 757 |
_supply[_node_id[s]] = k; |
807 | 758 |
_supply[_node_id[t]] = -k; |
808 | 759 |
return *this; |
809 | 760 |
} |
810 | 761 |
|
811 | 762 |
/// \brief Set the type of the supply constraints. |
812 | 763 |
/// |
813 | 764 |
/// This function sets the type of the supply/demand constraints. |
814 | 765 |
/// If it is not used before calling \ref run(), the \ref GEQ supply |
815 | 766 |
/// type will be used. |
816 | 767 |
/// |
817 | 768 |
/// For more information, see \ref SupplyType. |
818 | 769 |
/// |
819 | 770 |
/// \return <tt>(*this)</tt> |
820 | 771 |
NetworkSimplex& supplyType(SupplyType supply_type) { |
821 | 772 |
_stype = supply_type; |
822 | 773 |
return *this; |
823 | 774 |
} |
824 | 775 |
|
825 | 776 |
/// @} |
826 | 777 |
|
827 | 778 |
/// \name Execution Control |
828 | 779 |
/// The algorithm can be executed using \ref run(). |
829 | 780 |
|
830 | 781 |
/// @{ |
831 | 782 |
|
832 | 783 |
/// \brief Run the algorithm. |
833 | 784 |
/// |
834 | 785 |
/// This function runs the algorithm. |
835 | 786 |
/// The paramters can be specified using functions \ref lowerMap(), |
836 | 787 |
/// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply(), |
837 | 788 |
/// \ref supplyType(). |
838 | 789 |
/// For example, |
839 | 790 |
/// \code |
840 | 791 |
/// NetworkSimplex<ListDigraph> ns(graph); |
841 | 792 |
/// ns.lowerMap(lower).upperMap(upper).costMap(cost) |
842 | 793 |
/// .supplyMap(sup).run(); |
843 | 794 |
/// \endcode |
844 | 795 |
/// |
845 |
/// This function can be called more than once. All the parameters |
|
846 |
/// that have been given are kept for the next call, unless |
|
847 |
/// \ref reset() is called, thus only the modified parameters |
|
848 |
/// have to be set again. See \ref reset() for examples. |
|
849 |
/// However, the underlying digraph must not be modified after this |
|
850 |
/// class have been constructed, since it copies and extends the graph. |
|
796 |
/// This function can be called more than once. All the given parameters |
|
797 |
/// are kept for the next call, unless \ref resetParams() or \ref reset() |
|
798 |
/// is used, thus only the modified parameters have to be set again. |
|
799 |
/// If the underlying digraph was also modified after the construction |
|
800 |
/// of the class (or the last \ref reset() call), then the \ref reset() |
|
801 |
/// function must be called. |
|
851 | 802 |
/// |
852 | 803 |
/// \param pivot_rule The pivot rule that will be used during the |
853 | 804 |
/// algorithm. For more information, see \ref PivotRule. |
854 | 805 |
/// |
855 | 806 |
/// \return \c INFEASIBLE if no feasible flow exists, |
856 | 807 |
/// \n \c OPTIMAL if the problem has optimal solution |
857 | 808 |
/// (i.e. it is feasible and bounded), and the algorithm has found |
858 | 809 |
/// optimal flow and node potentials (primal and dual solutions), |
859 | 810 |
/// \n \c UNBOUNDED if the objective function of the problem is |
860 | 811 |
/// unbounded, i.e. there is a directed cycle having negative total |
861 | 812 |
/// cost and infinite upper bound. |
862 | 813 |
/// |
863 | 814 |
/// \see ProblemType, PivotRule |
815 |
/// \see resetParams(), reset() |
|
864 | 816 |
ProblemType run(PivotRule pivot_rule = BLOCK_SEARCH) { |
865 | 817 |
if (!init()) return INFEASIBLE; |
866 | 818 |
return start(pivot_rule); |
867 | 819 |
} |
868 | 820 |
|
869 | 821 |
/// \brief Reset all the parameters that have been given before. |
870 | 822 |
/// |
871 | 823 |
/// This function resets all the paramaters that have been given |
872 | 824 |
/// before using functions \ref lowerMap(), \ref upperMap(), |
873 | 825 |
/// \ref costMap(), \ref supplyMap(), \ref stSupply(), \ref supplyType(). |
874 | 826 |
/// |
875 |
/// It is useful for multiple run() calls. If this function is not |
|
876 |
/// used, all the parameters given before are kept for the next |
|
877 |
/// \ref run() call. |
|
878 |
/// However, the underlying digraph must not be modified after this |
|
879 |
/// |
|
827 |
/// It is useful for multiple \ref run() calls. Basically, all the given |
|
828 |
/// parameters are kept for the next \ref run() call, unless |
|
829 |
/// \ref resetParams() or \ref reset() is used. |
|
830 |
/// If the underlying digraph was also modified after the construction |
|
831 |
/// of the class or the last \ref reset() call, then the \ref reset() |
|
832 |
/// function must be used, otherwise \ref resetParams() is sufficient. |
|
880 | 833 |
/// |
881 | 834 |
/// For example, |
882 | 835 |
/// \code |
883 | 836 |
/// NetworkSimplex<ListDigraph> ns(graph); |
884 | 837 |
/// |
885 | 838 |
/// // First run |
886 | 839 |
/// ns.lowerMap(lower).upperMap(upper).costMap(cost) |
887 | 840 |
/// .supplyMap(sup).run(); |
888 | 841 |
/// |
889 |
/// // Run again with modified cost map ( |
|
842 |
/// // Run again with modified cost map (resetParams() is not called, |
|
890 | 843 |
/// // so only the cost map have to be set again) |
891 | 844 |
/// cost[e] += 100; |
892 | 845 |
/// ns.costMap(cost).run(); |
893 | 846 |
/// |
894 |
/// // Run again from scratch using |
|
847 |
/// // Run again from scratch using resetParams() |
|
895 | 848 |
/// // (the lower bounds will be set to zero on all arcs) |
896 |
/// ns. |
|
849 |
/// ns.resetParams(); |
|
897 | 850 |
/// ns.upperMap(capacity).costMap(cost) |
898 | 851 |
/// .supplyMap(sup).run(); |
899 | 852 |
/// \endcode |
900 | 853 |
/// |
901 | 854 |
/// \return <tt>(*this)</tt> |
902 |
|
|
855 |
/// |
|
856 |
/// \see reset(), run() |
|
857 |
NetworkSimplex& resetParams() { |
|
903 | 858 |
for (int i = 0; i != _node_num; ++i) { |
904 | 859 |
_supply[i] = 0; |
905 | 860 |
} |
906 | 861 |
for (int i = 0; i != _arc_num; ++i) { |
907 | 862 |
_lower[i] = 0; |
908 | 863 |
_upper[i] = INF; |
909 | 864 |
_cost[i] = 1; |
910 | 865 |
} |
911 | 866 |
_have_lower = false; |
912 | 867 |
_stype = GEQ; |
913 | 868 |
return *this; |
914 | 869 |
} |
915 | 870 |
|
871 |
/// \brief Reset the internal data structures and all the parameters |
|
872 |
/// that have been given before. |
|
873 |
/// |
|
874 |
/// This function resets the internal data structures and all the |
|
875 |
/// paramaters that have been given before using functions \ref lowerMap(), |
|
876 |
/// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply(), |
|
877 |
/// \ref supplyType(). |
|
878 |
/// |
|
879 |
/// It is useful for multiple \ref run() calls. Basically, all the given |
|
880 |
/// parameters are kept for the next \ref run() call, unless |
|
881 |
/// \ref resetParams() or \ref reset() is used. |
|
882 |
/// If the underlying digraph was also modified after the construction |
|
883 |
/// of the class or the last \ref reset() call, then the \ref reset() |
|
884 |
/// function must be used, otherwise \ref resetParams() is sufficient. |
|
885 |
/// |
|
886 |
/// See \ref resetParams() for examples. |
|
887 |
/// |
|
888 |
/// \return <tt>(*this)</tt> |
|
889 |
/// |
|
890 |
/// \see resetParams(), run() |
|
891 |
NetworkSimplex& reset() { |
|
892 |
// Resize vectors |
|
893 |
_node_num = countNodes(_graph); |
|
894 |
_arc_num = countArcs(_graph); |
|
895 |
int all_node_num = _node_num + 1; |
|
896 |
int max_arc_num = _arc_num + 2 * _node_num; |
|
897 |
|
|
898 |
_source.resize(max_arc_num); |
|
899 |
_target.resize(max_arc_num); |
|
900 |
|
|
901 |
_lower.resize(_arc_num); |
|
902 |
_upper.resize(_arc_num); |
|
903 |
_cap.resize(max_arc_num); |
|
904 |
_cost.resize(max_arc_num); |
|
905 |
_supply.resize(all_node_num); |
|
906 |
_flow.resize(max_arc_num); |
|
907 |
_pi.resize(all_node_num); |
|
908 |
|
|
909 |
_parent.resize(all_node_num); |
|
910 |
_pred.resize(all_node_num); |
|
911 |
_forward.resize(all_node_num); |
|
912 |
_thread.resize(all_node_num); |
|
913 |
_rev_thread.resize(all_node_num); |
|
914 |
_succ_num.resize(all_node_num); |
|
915 |
_last_succ.resize(all_node_num); |
|
916 |
_state.resize(max_arc_num); |
|
917 |
|
|
918 |
// Copy the graph |
|
919 |
int i = 0; |
|
920 |
for (NodeIt n(_graph); n != INVALID; ++n, ++i) { |
|
921 |
_node_id[n] = i; |
|
922 |
} |
|
923 |
if (_arc_mixing) { |
|
924 |
// Store the arcs in a mixed order |
|
925 |
int k = std::max(int(std::sqrt(double(_arc_num))), 10); |
|
926 |
int i = 0, j = 0; |
|
927 |
for (ArcIt a(_graph); a != INVALID; ++a) { |
|
928 |
_arc_id[a] = i; |
|
929 |
_source[i] = _node_id[_graph.source(a)]; |
|
930 |
_target[i] = _node_id[_graph.target(a)]; |
|
931 |
if ((i += k) >= _arc_num) i = ++j; |
|
932 |
} |
|
933 |
} else { |
|
934 |
// Store the arcs in the original order |
|
935 |
int i = 0; |
|
936 |
for (ArcIt a(_graph); a != INVALID; ++a, ++i) { |
|
937 |
_arc_id[a] = i; |
|
938 |
_source[i] = _node_id[_graph.source(a)]; |
|
939 |
_target[i] = _node_id[_graph.target(a)]; |
|
940 |
} |
|
941 |
} |
|
942 |
|
|
943 |
// Reset parameters |
|
944 |
resetParams(); |
|
945 |
return *this; |
|
946 |
} |
|
947 |
|
|
916 | 948 |
/// @} |
917 | 949 |
|
918 | 950 |
/// \name Query Functions |
919 | 951 |
/// The results of the algorithm can be obtained using these |
920 | 952 |
/// functions.\n |
921 | 953 |
/// The \ref run() function must be called before using them. |
922 | 954 |
|
923 | 955 |
/// @{ |
924 | 956 |
|
925 | 957 |
/// \brief Return the total cost of the found flow. |
926 | 958 |
/// |
927 | 959 |
/// This function returns the total cost of the found flow. |
928 | 960 |
/// Its complexity is O(e). |
929 | 961 |
/// |
930 | 962 |
/// \note The return type of the function can be specified as a |
931 | 963 |
/// template parameter. For example, |
932 | 964 |
/// \code |
933 | 965 |
/// ns.totalCost<double>(); |
934 | 966 |
/// \endcode |
935 | 967 |
/// It is useful if the total cost cannot be stored in the \c Cost |
936 | 968 |
/// type of the algorithm, which is the default return type of the |
937 | 969 |
/// function. |
938 | 970 |
/// |
939 | 971 |
/// \pre \ref run() must be called before using this function. |
940 | 972 |
template <typename Number> |
941 | 973 |
Number totalCost() const { |
942 | 974 |
Number c = 0; |
943 | 975 |
for (ArcIt a(_graph); a != INVALID; ++a) { |
944 | 976 |
int i = _arc_id[a]; |
945 | 977 |
c += Number(_flow[i]) * Number(_cost[i]); |
946 | 978 |
} |
947 | 979 |
return c; |
948 | 980 |
} |
949 | 981 |
|
950 | 982 |
#ifndef DOXYGEN |
951 | 983 |
Cost totalCost() const { |
952 | 984 |
return totalCost<Cost>(); |
953 | 985 |
} |
954 | 986 |
#endif |
955 | 987 |
|
956 | 988 |
/// \brief Return the flow on the given arc. |
957 | 989 |
/// |
958 | 990 |
/// This function returns the flow on the given arc. |
959 | 991 |
/// |
960 | 992 |
/// \pre \ref run() must be called before using this function. |
961 | 993 |
Value flow(const Arc& a) const { |
962 | 994 |
return _flow[_arc_id[a]]; |
963 | 995 |
} |
964 | 996 |
|
965 | 997 |
/// \brief Return the flow map (the primal solution). |
966 | 998 |
/// |
967 | 999 |
/// This function copies the flow value on each arc into the given |
968 | 1000 |
/// map. The \c Value type of the algorithm must be convertible to |
969 | 1001 |
/// the \c Value type of the map. |
970 | 1002 |
/// |
971 | 1003 |
/// \pre \ref run() must be called before using this function. |
972 | 1004 |
template <typename FlowMap> |
973 | 1005 |
void flowMap(FlowMap &map) const { |
974 | 1006 |
for (ArcIt a(_graph); a != INVALID; ++a) { |
975 | 1007 |
map.set(a, _flow[_arc_id[a]]); |
976 | 1008 |
} |
977 | 1009 |
} |
978 | 1010 |
|
979 | 1011 |
/// \brief Return the potential (dual value) of the given node. |
980 | 1012 |
/// |
981 | 1013 |
/// This function returns the potential (dual value) of the |
982 | 1014 |
/// given node. |
983 | 1015 |
/// |
984 | 1016 |
/// \pre \ref run() must be called before using this function. |
985 | 1017 |
Cost potential(const Node& n) const { |
986 | 1018 |
return _pi[_node_id[n]]; |
987 | 1019 |
} |
988 | 1020 |
|
989 | 1021 |
/// \brief Return the potential map (the dual solution). |
990 | 1022 |
/// |
991 | 1023 |
/// This function copies the potential (dual value) of each node |
992 | 1024 |
/// into the given map. |
993 | 1025 |
/// The \c Cost type of the algorithm must be convertible to the |
994 | 1026 |
/// \c Value type of the map. |
995 | 1027 |
/// |
996 | 1028 |
/// \pre \ref run() must be called before using this function. |
997 | 1029 |
template <typename PotentialMap> |
998 | 1030 |
void potentialMap(PotentialMap &map) const { |
999 | 1031 |
for (NodeIt n(_graph); n != INVALID; ++n) { |
1000 | 1032 |
map.set(n, _pi[_node_id[n]]); |
1001 | 1033 |
} |
1002 | 1034 |
} |
1003 | 1035 |
|
1004 | 1036 |
/// @} |
1005 | 1037 |
|
1006 | 1038 |
private: |
1007 | 1039 |
|
1008 | 1040 |
// Initialize internal data structures |
1009 | 1041 |
bool init() { |
1010 | 1042 |
if (_node_num == 0) return false; |
1011 | 1043 |
... | ... |
@@ -64,416 +64,416 @@ |
64 | 64 |
" 4 7 80 2 0 0 0\n" |
65 | 65 |
" 4 8 110 3 0 0 0\n" |
66 | 66 |
" 5 7 60 14 0 0 0\n" |
67 | 67 |
" 5 11 120 12 0 0 0\n" |
68 | 68 |
" 6 3 0 3 0 0 0\n" |
69 | 69 |
" 6 9 140 4 0 0 0\n" |
70 | 70 |
" 6 10 90 8 0 0 0\n" |
71 | 71 |
" 7 1 30 5 0 0 -5\n" |
72 | 72 |
" 8 12 60 16 0 4 3\n" |
73 | 73 |
" 9 12 50 6 0 0 0\n" |
74 | 74 |
"10 12 70 13 0 5 2\n" |
75 | 75 |
"10 2 100 7 0 0 0\n" |
76 | 76 |
"10 7 60 10 0 0 -3\n" |
77 | 77 |
"11 10 20 14 0 6 -20\n" |
78 | 78 |
"12 11 30 10 0 0 -10\n" |
79 | 79 |
"\n" |
80 | 80 |
"@attributes\n" |
81 | 81 |
"source 1\n" |
82 | 82 |
"target 12\n"; |
83 | 83 |
|
84 | 84 |
char test_neg1_lgf[] = |
85 | 85 |
"@nodes\n" |
86 | 86 |
"label sup\n" |
87 | 87 |
" 1 100\n" |
88 | 88 |
" 2 0\n" |
89 | 89 |
" 3 0\n" |
90 | 90 |
" 4 -100\n" |
91 | 91 |
" 5 0\n" |
92 | 92 |
" 6 0\n" |
93 | 93 |
" 7 0\n" |
94 | 94 |
"@arcs\n" |
95 | 95 |
" cost low1 low2\n" |
96 | 96 |
"1 2 100 0 0\n" |
97 | 97 |
"1 3 30 0 0\n" |
98 | 98 |
"2 4 20 0 0\n" |
99 | 99 |
"3 4 80 0 0\n" |
100 | 100 |
"3 2 50 0 0\n" |
101 | 101 |
"5 3 10 0 0\n" |
102 | 102 |
"5 6 80 0 1000\n" |
103 | 103 |
"6 7 30 0 -1000\n" |
104 | 104 |
"7 5 -120 0 0\n"; |
105 | 105 |
|
106 | 106 |
char test_neg2_lgf[] = |
107 | 107 |
"@nodes\n" |
108 | 108 |
"label sup\n" |
109 | 109 |
" 1 100\n" |
110 | 110 |
" 2 -300\n" |
111 | 111 |
"@arcs\n" |
112 | 112 |
" cost\n" |
113 | 113 |
"1 2 -1\n"; |
114 | 114 |
|
115 | 115 |
|
116 | 116 |
// Test data |
117 | 117 |
typedef ListDigraph Digraph; |
118 | 118 |
DIGRAPH_TYPEDEFS(ListDigraph); |
119 | 119 |
|
120 | 120 |
Digraph gr; |
121 | 121 |
Digraph::ArcMap<int> c(gr), l1(gr), l2(gr), l3(gr), u(gr); |
122 | 122 |
Digraph::NodeMap<int> s1(gr), s2(gr), s3(gr), s4(gr), s5(gr), s6(gr); |
123 | 123 |
ConstMap<Arc, int> cc(1), cu(std::numeric_limits<int>::max()); |
124 | 124 |
Node v, w; |
125 | 125 |
|
126 | 126 |
Digraph neg1_gr; |
127 | 127 |
Digraph::ArcMap<int> neg1_c(neg1_gr), neg1_l1(neg1_gr), neg1_l2(neg1_gr); |
128 | 128 |
ConstMap<Arc, int> neg1_u1(std::numeric_limits<int>::max()), neg1_u2(5000); |
129 | 129 |
Digraph::NodeMap<int> neg1_s(neg1_gr); |
130 | 130 |
|
131 | 131 |
Digraph neg2_gr; |
132 | 132 |
Digraph::ArcMap<int> neg2_c(neg2_gr); |
133 | 133 |
ConstMap<Arc, int> neg2_l(0), neg2_u(1000); |
134 | 134 |
Digraph::NodeMap<int> neg2_s(neg2_gr); |
135 | 135 |
|
136 | 136 |
|
137 | 137 |
enum SupplyType { |
138 | 138 |
EQ, |
139 | 139 |
GEQ, |
140 | 140 |
LEQ |
141 | 141 |
}; |
142 | 142 |
|
143 | 143 |
|
144 | 144 |
// Check the interface of an MCF algorithm |
145 | 145 |
template <typename GR, typename Value, typename Cost> |
146 | 146 |
class McfClassConcept |
147 | 147 |
{ |
148 | 148 |
public: |
149 | 149 |
|
150 | 150 |
template <typename MCF> |
151 | 151 |
struct Constraints { |
152 | 152 |
void constraints() { |
153 | 153 |
checkConcept<concepts::Digraph, GR>(); |
154 | 154 |
|
155 | 155 |
const Constraints& me = *this; |
156 | 156 |
|
157 | 157 |
MCF mcf(me.g); |
158 | 158 |
const MCF& const_mcf = mcf; |
159 | 159 |
|
160 |
b = mcf.reset() |
|
160 |
b = mcf.reset().resetParams() |
|
161 | 161 |
.lowerMap(me.lower) |
162 | 162 |
.upperMap(me.upper) |
163 | 163 |
.costMap(me.cost) |
164 | 164 |
.supplyMap(me.sup) |
165 | 165 |
.stSupply(me.n, me.n, me.k) |
166 | 166 |
.run(); |
167 | 167 |
|
168 | 168 |
c = const_mcf.totalCost(); |
169 | 169 |
x = const_mcf.template totalCost<double>(); |
170 | 170 |
v = const_mcf.flow(me.a); |
171 | 171 |
c = const_mcf.potential(me.n); |
172 | 172 |
const_mcf.flowMap(fm); |
173 | 173 |
const_mcf.potentialMap(pm); |
174 | 174 |
} |
175 | 175 |
|
176 | 176 |
typedef typename GR::Node Node; |
177 | 177 |
typedef typename GR::Arc Arc; |
178 | 178 |
typedef concepts::ReadMap<Node, Value> NM; |
179 | 179 |
typedef concepts::ReadMap<Arc, Value> VAM; |
180 | 180 |
typedef concepts::ReadMap<Arc, Cost> CAM; |
181 | 181 |
typedef concepts::WriteMap<Arc, Value> FlowMap; |
182 | 182 |
typedef concepts::WriteMap<Node, Cost> PotMap; |
183 | 183 |
|
184 | 184 |
GR g; |
185 | 185 |
VAM lower; |
186 | 186 |
VAM upper; |
187 | 187 |
CAM cost; |
188 | 188 |
NM sup; |
189 | 189 |
Node n; |
190 | 190 |
Arc a; |
191 | 191 |
Value k; |
192 | 192 |
|
193 | 193 |
FlowMap fm; |
194 | 194 |
PotMap pm; |
195 | 195 |
bool b; |
196 | 196 |
double x; |
197 | 197 |
typename MCF::Value v; |
198 | 198 |
typename MCF::Cost c; |
199 | 199 |
}; |
200 | 200 |
|
201 | 201 |
}; |
202 | 202 |
|
203 | 203 |
|
204 | 204 |
// Check the feasibility of the given flow (primal soluiton) |
205 | 205 |
template < typename GR, typename LM, typename UM, |
206 | 206 |
typename SM, typename FM > |
207 | 207 |
bool checkFlow( const GR& gr, const LM& lower, const UM& upper, |
208 | 208 |
const SM& supply, const FM& flow, |
209 | 209 |
SupplyType type = EQ ) |
210 | 210 |
{ |
211 | 211 |
TEMPLATE_DIGRAPH_TYPEDEFS(GR); |
212 | 212 |
|
213 | 213 |
for (ArcIt e(gr); e != INVALID; ++e) { |
214 | 214 |
if (flow[e] < lower[e] || flow[e] > upper[e]) return false; |
215 | 215 |
} |
216 | 216 |
|
217 | 217 |
for (NodeIt n(gr); n != INVALID; ++n) { |
218 | 218 |
typename SM::Value sum = 0; |
219 | 219 |
for (OutArcIt e(gr, n); e != INVALID; ++e) |
220 | 220 |
sum += flow[e]; |
221 | 221 |
for (InArcIt e(gr, n); e != INVALID; ++e) |
222 | 222 |
sum -= flow[e]; |
223 | 223 |
bool b = (type == EQ && sum == supply[n]) || |
224 | 224 |
(type == GEQ && sum >= supply[n]) || |
225 | 225 |
(type == LEQ && sum <= supply[n]); |
226 | 226 |
if (!b) return false; |
227 | 227 |
} |
228 | 228 |
|
229 | 229 |
return true; |
230 | 230 |
} |
231 | 231 |
|
232 | 232 |
// Check the feasibility of the given potentials (dual soluiton) |
233 | 233 |
// using the "Complementary Slackness" optimality condition |
234 | 234 |
template < typename GR, typename LM, typename UM, |
235 | 235 |
typename CM, typename SM, typename FM, typename PM > |
236 | 236 |
bool checkPotential( const GR& gr, const LM& lower, const UM& upper, |
237 | 237 |
const CM& cost, const SM& supply, const FM& flow, |
238 | 238 |
const PM& pi, SupplyType type ) |
239 | 239 |
{ |
240 | 240 |
TEMPLATE_DIGRAPH_TYPEDEFS(GR); |
241 | 241 |
|
242 | 242 |
bool opt = true; |
243 | 243 |
for (ArcIt e(gr); opt && e != INVALID; ++e) { |
244 | 244 |
typename CM::Value red_cost = |
245 | 245 |
cost[e] + pi[gr.source(e)] - pi[gr.target(e)]; |
246 | 246 |
opt = red_cost == 0 || |
247 | 247 |
(red_cost > 0 && flow[e] == lower[e]) || |
248 | 248 |
(red_cost < 0 && flow[e] == upper[e]); |
249 | 249 |
} |
250 | 250 |
|
251 | 251 |
for (NodeIt n(gr); opt && n != INVALID; ++n) { |
252 | 252 |
typename SM::Value sum = 0; |
253 | 253 |
for (OutArcIt e(gr, n); e != INVALID; ++e) |
254 | 254 |
sum += flow[e]; |
255 | 255 |
for (InArcIt e(gr, n); e != INVALID; ++e) |
256 | 256 |
sum -= flow[e]; |
257 | 257 |
if (type != LEQ) { |
258 | 258 |
opt = (pi[n] <= 0) && (sum == supply[n] || pi[n] == 0); |
259 | 259 |
} else { |
260 | 260 |
opt = (pi[n] >= 0) && (sum == supply[n] || pi[n] == 0); |
261 | 261 |
} |
262 | 262 |
} |
263 | 263 |
|
264 | 264 |
return opt; |
265 | 265 |
} |
266 | 266 |
|
267 | 267 |
// Check whether the dual cost is equal to the primal cost |
268 | 268 |
template < typename GR, typename LM, typename UM, |
269 | 269 |
typename CM, typename SM, typename PM > |
270 | 270 |
bool checkDualCost( const GR& gr, const LM& lower, const UM& upper, |
271 | 271 |
const CM& cost, const SM& supply, const PM& pi, |
272 | 272 |
typename CM::Value total ) |
273 | 273 |
{ |
274 | 274 |
TEMPLATE_DIGRAPH_TYPEDEFS(GR); |
275 | 275 |
|
276 | 276 |
typename CM::Value dual_cost = 0; |
277 | 277 |
SM red_supply(gr); |
278 | 278 |
for (NodeIt n(gr); n != INVALID; ++n) { |
279 | 279 |
red_supply[n] = supply[n]; |
280 | 280 |
} |
281 | 281 |
for (ArcIt a(gr); a != INVALID; ++a) { |
282 | 282 |
if (lower[a] != 0) { |
283 | 283 |
dual_cost += lower[a] * cost[a]; |
284 | 284 |
red_supply[gr.source(a)] -= lower[a]; |
285 | 285 |
red_supply[gr.target(a)] += lower[a]; |
286 | 286 |
} |
287 | 287 |
} |
288 | 288 |
|
289 | 289 |
for (NodeIt n(gr); n != INVALID; ++n) { |
290 | 290 |
dual_cost -= red_supply[n] * pi[n]; |
291 | 291 |
} |
292 | 292 |
for (ArcIt a(gr); a != INVALID; ++a) { |
293 | 293 |
typename CM::Value red_cost = |
294 | 294 |
cost[a] + pi[gr.source(a)] - pi[gr.target(a)]; |
295 | 295 |
dual_cost -= (upper[a] - lower[a]) * std::max(-red_cost, 0); |
296 | 296 |
} |
297 | 297 |
|
298 | 298 |
return dual_cost == total; |
299 | 299 |
} |
300 | 300 |
|
301 | 301 |
// Run a minimum cost flow algorithm and check the results |
302 | 302 |
template < typename MCF, typename GR, |
303 | 303 |
typename LM, typename UM, |
304 | 304 |
typename CM, typename SM, |
305 | 305 |
typename PT > |
306 | 306 |
void checkMcf( const MCF& mcf, PT mcf_result, |
307 | 307 |
const GR& gr, const LM& lower, const UM& upper, |
308 | 308 |
const CM& cost, const SM& supply, |
309 | 309 |
PT result, bool optimal, typename CM::Value total, |
310 | 310 |
const std::string &test_id = "", |
311 | 311 |
SupplyType type = EQ ) |
312 | 312 |
{ |
313 | 313 |
check(mcf_result == result, "Wrong result " + test_id); |
314 | 314 |
if (optimal) { |
315 | 315 |
typename GR::template ArcMap<typename SM::Value> flow(gr); |
316 | 316 |
typename GR::template NodeMap<typename CM::Value> pi(gr); |
317 | 317 |
mcf.flowMap(flow); |
318 | 318 |
mcf.potentialMap(pi); |
319 | 319 |
check(checkFlow(gr, lower, upper, supply, flow, type), |
320 | 320 |
"The flow is not feasible " + test_id); |
321 | 321 |
check(mcf.totalCost() == total, "The flow is not optimal " + test_id); |
322 | 322 |
check(checkPotential(gr, lower, upper, cost, supply, flow, pi, type), |
323 | 323 |
"Wrong potentials " + test_id); |
324 | 324 |
check(checkDualCost(gr, lower, upper, cost, supply, pi, total), |
325 | 325 |
"Wrong dual cost " + test_id); |
326 | 326 |
} |
327 | 327 |
} |
328 | 328 |
|
329 | 329 |
template < typename MCF, typename Param > |
330 | 330 |
void runMcfGeqTests( Param param, |
331 | 331 |
const std::string &test_str = "", |
332 | 332 |
bool full_neg_cost_support = false ) |
333 | 333 |
{ |
334 | 334 |
MCF mcf1(gr), mcf2(neg1_gr), mcf3(neg2_gr); |
335 | 335 |
|
336 | 336 |
// Basic tests |
337 | 337 |
mcf1.upperMap(u).costMap(c).supplyMap(s1); |
338 | 338 |
checkMcf(mcf1, mcf1.run(param), gr, l1, u, c, s1, |
339 | 339 |
mcf1.OPTIMAL, true, 5240, test_str + "-1"); |
340 | 340 |
mcf1.stSupply(v, w, 27); |
341 | 341 |
checkMcf(mcf1, mcf1.run(param), gr, l1, u, c, s2, |
342 | 342 |
mcf1.OPTIMAL, true, 7620, test_str + "-2"); |
343 | 343 |
mcf1.lowerMap(l2).supplyMap(s1); |
344 | 344 |
checkMcf(mcf1, mcf1.run(param), gr, l2, u, c, s1, |
345 | 345 |
mcf1.OPTIMAL, true, 5970, test_str + "-3"); |
346 | 346 |
mcf1.stSupply(v, w, 27); |
347 | 347 |
checkMcf(mcf1, mcf1.run(param), gr, l2, u, c, s2, |
348 | 348 |
mcf1.OPTIMAL, true, 8010, test_str + "-4"); |
349 |
mcf1. |
|
349 |
mcf1.resetParams().supplyMap(s1); |
|
350 | 350 |
checkMcf(mcf1, mcf1.run(param), gr, l1, cu, cc, s1, |
351 | 351 |
mcf1.OPTIMAL, true, 74, test_str + "-5"); |
352 | 352 |
mcf1.lowerMap(l2).stSupply(v, w, 27); |
353 | 353 |
checkMcf(mcf1, mcf1.run(param), gr, l2, cu, cc, s2, |
354 | 354 |
mcf1.OPTIMAL, true, 94, test_str + "-6"); |
355 | 355 |
mcf1.reset(); |
356 | 356 |
checkMcf(mcf1, mcf1.run(param), gr, l1, cu, cc, s3, |
357 | 357 |
mcf1.OPTIMAL, true, 0, test_str + "-7"); |
358 | 358 |
mcf1.lowerMap(l2).upperMap(u); |
359 | 359 |
checkMcf(mcf1, mcf1.run(param), gr, l2, u, cc, s3, |
360 | 360 |
mcf1.INFEASIBLE, false, 0, test_str + "-8"); |
361 | 361 |
mcf1.lowerMap(l3).upperMap(u).costMap(c).supplyMap(s4); |
362 | 362 |
checkMcf(mcf1, mcf1.run(param), gr, l3, u, c, s4, |
363 | 363 |
mcf1.OPTIMAL, true, 6360, test_str + "-9"); |
364 | 364 |
|
365 | 365 |
// Tests for the GEQ form |
366 |
mcf1. |
|
366 |
mcf1.resetParams().upperMap(u).costMap(c).supplyMap(s5); |
|
367 | 367 |
checkMcf(mcf1, mcf1.run(param), gr, l1, u, c, s5, |
368 | 368 |
mcf1.OPTIMAL, true, 3530, test_str + "-10", GEQ); |
369 | 369 |
mcf1.lowerMap(l2); |
370 | 370 |
checkMcf(mcf1, mcf1.run(param), gr, l2, u, c, s5, |
371 | 371 |
mcf1.OPTIMAL, true, 4540, test_str + "-11", GEQ); |
372 | 372 |
mcf1.supplyMap(s6); |
373 | 373 |
checkMcf(mcf1, mcf1.run(param), gr, l2, u, c, s6, |
374 | 374 |
mcf1.INFEASIBLE, false, 0, test_str + "-12", GEQ); |
375 | 375 |
|
376 | 376 |
// Tests with negative costs |
377 | 377 |
mcf2.lowerMap(neg1_l1).costMap(neg1_c).supplyMap(neg1_s); |
378 | 378 |
checkMcf(mcf2, mcf2.run(param), neg1_gr, neg1_l1, neg1_u1, neg1_c, neg1_s, |
379 | 379 |
mcf2.UNBOUNDED, false, 0, test_str + "-13"); |
380 | 380 |
mcf2.upperMap(neg1_u2); |
381 | 381 |
checkMcf(mcf2, mcf2.run(param), neg1_gr, neg1_l1, neg1_u2, neg1_c, neg1_s, |
382 | 382 |
mcf2.OPTIMAL, true, -40000, test_str + "-14"); |
383 |
mcf2. |
|
383 |
mcf2.resetParams().lowerMap(neg1_l2).costMap(neg1_c).supplyMap(neg1_s); |
|
384 | 384 |
checkMcf(mcf2, mcf2.run(param), neg1_gr, neg1_l2, neg1_u1, neg1_c, neg1_s, |
385 | 385 |
mcf2.UNBOUNDED, false, 0, test_str + "-15"); |
386 | 386 |
|
387 | 387 |
mcf3.costMap(neg2_c).supplyMap(neg2_s); |
388 | 388 |
if (full_neg_cost_support) { |
389 | 389 |
checkMcf(mcf3, mcf3.run(param), neg2_gr, neg2_l, neg2_u, neg2_c, neg2_s, |
390 | 390 |
mcf3.OPTIMAL, true, -300, test_str + "-16", GEQ); |
391 | 391 |
} else { |
392 | 392 |
checkMcf(mcf3, mcf3.run(param), neg2_gr, neg2_l, neg2_u, neg2_c, neg2_s, |
393 | 393 |
mcf3.UNBOUNDED, false, 0, test_str + "-17", GEQ); |
394 | 394 |
} |
395 | 395 |
mcf3.upperMap(neg2_u); |
396 | 396 |
checkMcf(mcf3, mcf3.run(param), neg2_gr, neg2_l, neg2_u, neg2_c, neg2_s, |
397 | 397 |
mcf3.OPTIMAL, true, -300, test_str + "-18", GEQ); |
398 | 398 |
} |
399 | 399 |
|
400 | 400 |
template < typename MCF, typename Param > |
401 | 401 |
void runMcfLeqTests( Param param, |
402 | 402 |
const std::string &test_str = "" ) |
403 | 403 |
{ |
404 | 404 |
// Tests for the LEQ form |
405 | 405 |
MCF mcf1(gr); |
406 | 406 |
mcf1.supplyType(mcf1.LEQ); |
407 | 407 |
mcf1.upperMap(u).costMap(c).supplyMap(s6); |
408 | 408 |
checkMcf(mcf1, mcf1.run(param), gr, l1, u, c, s6, |
409 | 409 |
mcf1.OPTIMAL, true, 5080, test_str + "-19", LEQ); |
410 | 410 |
mcf1.lowerMap(l2); |
411 | 411 |
checkMcf(mcf1, mcf1.run(param), gr, l2, u, c, s6, |
412 | 412 |
mcf1.OPTIMAL, true, 5930, test_str + "-20", LEQ); |
413 | 413 |
mcf1.supplyMap(s5); |
414 | 414 |
checkMcf(mcf1, mcf1.run(param), gr, l2, u, c, s5, |
415 | 415 |
mcf1.INFEASIBLE, false, 0, test_str + "-21", LEQ); |
416 | 416 |
} |
417 | 417 |
|
418 | 418 |
|
419 | 419 |
int main() |
420 | 420 |
{ |
421 | 421 |
// Read the test networks |
422 | 422 |
std::istringstream input(test_lgf); |
423 | 423 |
DigraphReader<Digraph>(gr, input) |
424 | 424 |
.arcMap("cost", c) |
425 | 425 |
.arcMap("cap", u) |
426 | 426 |
.arcMap("low1", l1) |
427 | 427 |
.arcMap("low2", l2) |
428 | 428 |
.arcMap("low3", l3) |
429 | 429 |
.nodeMap("sup1", s1) |
430 | 430 |
.nodeMap("sup2", s2) |
431 | 431 |
.nodeMap("sup3", s3) |
432 | 432 |
.nodeMap("sup4", s4) |
433 | 433 |
.nodeMap("sup5", s5) |
434 | 434 |
.nodeMap("sup6", s6) |
435 | 435 |
.node("source", v) |
436 | 436 |
.node("target", w) |
437 | 437 |
.run(); |
438 | 438 |
|
439 | 439 |
std::istringstream neg_inp1(test_neg1_lgf); |
440 | 440 |
DigraphReader<Digraph>(neg1_gr, neg_inp1) |
441 | 441 |
.arcMap("cost", neg1_c) |
442 | 442 |
.arcMap("low1", neg1_l1) |
443 | 443 |
.arcMap("low2", neg1_l2) |
444 | 444 |
.nodeMap("sup", neg1_s) |
445 | 445 |
.run(); |
446 | 446 |
|
447 | 447 |
std::istringstream neg_inp2(test_neg2_lgf); |
448 | 448 |
DigraphReader<Digraph>(neg2_gr, neg_inp2) |
449 | 449 |
.arcMap("cost", neg2_c) |
450 | 450 |
.nodeMap("sup", neg2_s) |
451 | 451 |
.run(); |
452 | 452 |
|
453 | 453 |
// Check the interface of NetworkSimplex |
454 | 454 |
{ |
455 | 455 |
typedef concepts::Digraph GR; |
456 | 456 |
checkConcept< McfClassConcept<GR, int, int>, |
457 | 457 |
NetworkSimplex<GR> >(); |
458 | 458 |
checkConcept< McfClassConcept<GR, double, double>, |
459 | 459 |
NetworkSimplex<GR, double> >(); |
460 | 460 |
checkConcept< McfClassConcept<GR, int, double>, |
461 | 461 |
NetworkSimplex<GR, int, double> >(); |
462 | 462 |
} |
463 | 463 |
|
464 | 464 |
// Check the interface of CapacityScaling |
465 | 465 |
{ |
466 | 466 |
typedef concepts::Digraph GR; |
467 | 467 |
checkConcept< McfClassConcept<GR, int, int>, |
468 | 468 |
CapacityScaling<GR> >(); |
469 | 469 |
checkConcept< McfClassConcept<GR, double, double>, |
470 | 470 |
CapacityScaling<GR, double> >(); |
471 | 471 |
checkConcept< McfClassConcept<GR, int, double>, |
472 | 472 |
CapacityScaling<GR, int, double> >(); |
473 | 473 |
typedef CapacityScaling<GR>:: |
474 | 474 |
SetHeap<concepts::Heap<int, RangeMap<int> > >::Create CAS; |
475 | 475 |
checkConcept< McfClassConcept<GR, int, int>, CAS >(); |
476 | 476 |
} |
477 | 477 |
|
478 | 478 |
// Check the interface of CostScaling |
479 | 479 |
{ |
0 comments (0 inline)