360 |
360 |
361 // Find next entering arc |
361 // Find next entering arc |
362 bool findEnteringArc() { |
362 bool findEnteringArc() { |
363 Cost c, min = 0; |
363 Cost c, min = 0; |
364 int cnt = _block_size; |
364 int cnt = _block_size; |
365 int e, min_arc = _next_arc; |
365 int e; |
366 for (e = _next_arc; e < _search_arc_num; ++e) { |
366 for (e = _next_arc; e < _search_arc_num; ++e) { |
367 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
367 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
368 if (c < min) { |
368 if (c < min) { |
369 min = c; |
369 min = c; |
370 min_arc = e; |
370 _in_arc = e; |
371 } |
371 } |
372 if (--cnt == 0) { |
372 if (--cnt == 0) { |
373 if (min < 0) break; |
373 if (min < 0) goto search_end; |
374 cnt = _block_size; |
374 cnt = _block_size; |
375 } |
375 } |
376 } |
376 } |
377 if (min == 0 || cnt > 0) { |
377 for (e = 0; e < _next_arc; ++e) { |
378 for (e = 0; e < _next_arc; ++e) { |
378 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
379 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
379 if (c < min) { |
380 if (c < min) { |
380 min = c; |
381 min = c; |
381 _in_arc = e; |
382 min_arc = e; |
382 } |
383 } |
383 if (--cnt == 0) { |
384 if (--cnt == 0) { |
384 if (min < 0) goto search_end; |
385 if (min < 0) break; |
385 cnt = _block_size; |
386 cnt = _block_size; |
|
387 } |
|
388 } |
386 } |
389 } |
387 } |
390 if (min >= 0) return false; |
388 if (min >= 0) return false; |
391 _in_arc = min_arc; |
389 |
|
390 search_end: |
392 _next_arc = e; |
391 _next_arc = e; |
393 return true; |
392 return true; |
394 } |
393 } |
395 |
394 |
396 }; //class BlockSearchPivotRule |
395 }; //class BlockSearchPivotRule |
424 _cost(ns._cost), _state(ns._state), _pi(ns._pi), |
423 _cost(ns._cost), _state(ns._state), _pi(ns._pi), |
425 _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num), |
424 _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num), |
426 _next_arc(0) |
425 _next_arc(0) |
427 { |
426 { |
428 // The main parameters of the pivot rule |
427 // The main parameters of the pivot rule |
429 const double LIST_LENGTH_FACTOR = 1.0; |
428 const double LIST_LENGTH_FACTOR = 0.25; |
430 const int MIN_LIST_LENGTH = 10; |
429 const int MIN_LIST_LENGTH = 10; |
431 const double MINOR_LIMIT_FACTOR = 0.1; |
430 const double MINOR_LIMIT_FACTOR = 0.1; |
432 const int MIN_MINOR_LIMIT = 3; |
431 const int MIN_MINOR_LIMIT = 3; |
433 |
432 |
434 _list_length = std::max( int(LIST_LENGTH_FACTOR * |
433 _list_length = std::max( int(LIST_LENGTH_FACTOR * |
441 } |
440 } |
442 |
441 |
443 /// Find next entering arc |
442 /// Find next entering arc |
444 bool findEnteringArc() { |
443 bool findEnteringArc() { |
445 Cost min, c; |
444 Cost min, c; |
446 int e, min_arc = _next_arc; |
445 int e; |
447 if (_curr_length > 0 && _minor_count < _minor_limit) { |
446 if (_curr_length > 0 && _minor_count < _minor_limit) { |
448 // Minor iteration: select the best eligible arc from the |
447 // Minor iteration: select the best eligible arc from the |
449 // current candidate list |
448 // current candidate list |
450 ++_minor_count; |
449 ++_minor_count; |
451 min = 0; |
450 min = 0; |
452 for (int i = 0; i < _curr_length; ++i) { |
451 for (int i = 0; i < _curr_length; ++i) { |
453 e = _candidates[i]; |
452 e = _candidates[i]; |
454 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
453 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
455 if (c < min) { |
454 if (c < min) { |
456 min = c; |
455 min = c; |
457 min_arc = e; |
456 _in_arc = e; |
458 } |
457 } |
459 if (c >= 0) { |
458 else if (c >= 0) { |
460 _candidates[i--] = _candidates[--_curr_length]; |
459 _candidates[i--] = _candidates[--_curr_length]; |
461 } |
460 } |
462 } |
461 } |
463 if (min < 0) { |
462 if (min < 0) return true; |
464 _in_arc = min_arc; |
|
465 return true; |
|
466 } |
|
467 } |
463 } |
468 |
464 |
469 // Major iteration: build a new candidate list |
465 // Major iteration: build a new candidate list |
470 min = 0; |
466 min = 0; |
471 _curr_length = 0; |
467 _curr_length = 0; |
473 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
469 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
474 if (c < 0) { |
470 if (c < 0) { |
475 _candidates[_curr_length++] = e; |
471 _candidates[_curr_length++] = e; |
476 if (c < min) { |
472 if (c < min) { |
477 min = c; |
473 min = c; |
478 min_arc = e; |
474 _in_arc = e; |
479 } |
475 } |
480 if (_curr_length == _list_length) break; |
476 if (_curr_length == _list_length) goto search_end; |
481 } |
477 } |
482 } |
478 } |
483 if (_curr_length < _list_length) { |
479 for (e = 0; e < _next_arc; ++e) { |
484 for (e = 0; e < _next_arc; ++e) { |
480 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
485 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
481 if (c < 0) { |
486 if (c < 0) { |
482 _candidates[_curr_length++] = e; |
487 _candidates[_curr_length++] = e; |
483 if (c < min) { |
488 if (c < min) { |
484 min = c; |
489 min = c; |
485 _in_arc = e; |
490 min_arc = e; |
|
491 } |
|
492 if (_curr_length == _list_length) break; |
|
493 } |
486 } |
|
487 if (_curr_length == _list_length) goto search_end; |
494 } |
488 } |
495 } |
489 } |
496 if (_curr_length == 0) return false; |
490 if (_curr_length == 0) return false; |
|
491 |
|
492 search_end: |
497 _minor_count = 1; |
493 _minor_count = 1; |
498 _in_arc = min_arc; |
|
499 _next_arc = e; |
494 _next_arc = e; |
500 return true; |
495 return true; |
501 } |
496 } |
502 |
497 |
503 }; //class CandidateListPivotRule |
498 }; //class CandidateListPivotRule |
545 _cost(ns._cost), _state(ns._state), _pi(ns._pi), |
540 _cost(ns._cost), _state(ns._state), _pi(ns._pi), |
546 _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num), |
541 _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num), |
547 _next_arc(0), _cand_cost(ns._search_arc_num), _sort_func(_cand_cost) |
542 _next_arc(0), _cand_cost(ns._search_arc_num), _sort_func(_cand_cost) |
548 { |
543 { |
549 // The main parameters of the pivot rule |
544 // The main parameters of the pivot rule |
550 const double BLOCK_SIZE_FACTOR = 1.5; |
545 const double BLOCK_SIZE_FACTOR = 1.0; |
551 const int MIN_BLOCK_SIZE = 10; |
546 const int MIN_BLOCK_SIZE = 10; |
552 const double HEAD_LENGTH_FACTOR = 0.1; |
547 const double HEAD_LENGTH_FACTOR = 0.1; |
553 const int MIN_HEAD_LENGTH = 3; |
548 const int MIN_HEAD_LENGTH = 3; |
554 |
549 |
555 _block_size = std::max( int(BLOCK_SIZE_FACTOR * |
550 _block_size = std::max( int(BLOCK_SIZE_FACTOR * |
574 } |
569 } |
575 } |
570 } |
576 |
571 |
577 // Extend the list |
572 // Extend the list |
578 int cnt = _block_size; |
573 int cnt = _block_size; |
579 int last_arc = 0; |
|
580 int limit = _head_length; |
574 int limit = _head_length; |
581 |
575 |
582 for (int e = _next_arc; e < _search_arc_num; ++e) { |
576 for (e = _next_arc; e < _search_arc_num; ++e) { |
583 _cand_cost[e] = _state[e] * |
577 _cand_cost[e] = _state[e] * |
584 (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
578 (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
585 if (_cand_cost[e] < 0) { |
579 if (_cand_cost[e] < 0) { |
586 _candidates[_curr_length++] = e; |
580 _candidates[_curr_length++] = e; |
587 last_arc = e; |
|
588 } |
581 } |
589 if (--cnt == 0) { |
582 if (--cnt == 0) { |
590 if (_curr_length > limit) break; |
583 if (_curr_length > limit) goto search_end; |
591 limit = 0; |
584 limit = 0; |
592 cnt = _block_size; |
585 cnt = _block_size; |
593 } |
586 } |
594 } |
587 } |
595 if (_curr_length <= limit) { |
588 for (e = 0; e < _next_arc; ++e) { |
596 for (int e = 0; e < _next_arc; ++e) { |
589 _cand_cost[e] = _state[e] * |
597 _cand_cost[e] = _state[e] * |
590 (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
598 (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
591 if (_cand_cost[e] < 0) { |
599 if (_cand_cost[e] < 0) { |
592 _candidates[_curr_length++] = e; |
600 _candidates[_curr_length++] = e; |
593 } |
601 last_arc = e; |
594 if (--cnt == 0) { |
602 } |
595 if (_curr_length > limit) goto search_end; |
603 if (--cnt == 0) { |
596 limit = 0; |
604 if (_curr_length > limit) break; |
597 cnt = _block_size; |
605 limit = 0; |
|
606 cnt = _block_size; |
|
607 } |
|
608 } |
598 } |
609 } |
599 } |
610 if (_curr_length == 0) return false; |
600 if (_curr_length == 0) return false; |
611 _next_arc = last_arc + 1; |
601 |
|
602 search_end: |
612 |
603 |
613 // Make heap of the candidate list (approximating a partial sort) |
604 // Make heap of the candidate list (approximating a partial sort) |
614 make_heap( _candidates.begin(), _candidates.begin() + _curr_length, |
605 make_heap( _candidates.begin(), _candidates.begin() + _curr_length, |
615 _sort_func ); |
606 _sort_func ); |
616 |
607 |
617 // Pop the first element of the heap |
608 // Pop the first element of the heap |
618 _in_arc = _candidates[0]; |
609 _in_arc = _candidates[0]; |
|
610 _next_arc = e; |
619 pop_heap( _candidates.begin(), _candidates.begin() + _curr_length, |
611 pop_heap( _candidates.begin(), _candidates.begin() + _curr_length, |
620 _sort_func ); |
612 _sort_func ); |
621 _curr_length = std::min(_head_length, _curr_length - 1); |
613 _curr_length = std::min(_head_length, _curr_length - 1); |
622 return true; |
614 return true; |
623 } |
615 } |
629 /// \brief Constructor. |
621 /// \brief Constructor. |
630 /// |
622 /// |
631 /// The constructor of the class. |
623 /// The constructor of the class. |
632 /// |
624 /// |
633 /// \param graph The digraph the algorithm runs on. |
625 /// \param graph The digraph the algorithm runs on. |
634 NetworkSimplex(const GR& graph) : |
626 /// \param arc_mixing Indicate if the arcs have to be stored in a |
|
627 /// mixed order in the internal data structure. |
|
628 /// In special cases, it could lead to better overall performance, |
|
629 /// but it is usually slower. Therefore it is disabled by default. |
|
630 NetworkSimplex(const GR& graph, bool arc_mixing = false) : |
635 _graph(graph), _node_id(graph), _arc_id(graph), |
631 _graph(graph), _node_id(graph), _arc_id(graph), |
636 INF(std::numeric_limits<Value>::has_infinity ? |
632 INF(std::numeric_limits<Value>::has_infinity ? |
637 std::numeric_limits<Value>::infinity() : |
633 std::numeric_limits<Value>::infinity() : |
638 std::numeric_limits<Value>::max()) |
634 std::numeric_limits<Value>::max()) |
639 { |
635 { |
667 _rev_thread.resize(all_node_num); |
663 _rev_thread.resize(all_node_num); |
668 _succ_num.resize(all_node_num); |
664 _succ_num.resize(all_node_num); |
669 _last_succ.resize(all_node_num); |
665 _last_succ.resize(all_node_num); |
670 _state.resize(max_arc_num); |
666 _state.resize(max_arc_num); |
671 |
667 |
672 // Copy the graph (store the arcs in a mixed order) |
668 // Copy the graph |
673 int i = 0; |
669 int i = 0; |
674 for (NodeIt n(_graph); n != INVALID; ++n, ++i) { |
670 for (NodeIt n(_graph); n != INVALID; ++n, ++i) { |
675 _node_id[n] = i; |
671 _node_id[n] = i; |
676 } |
672 } |
677 int k = std::max(int(std::sqrt(double(_arc_num))), 10); |
673 if (arc_mixing) { |
678 i = 0; |
674 // Store the arcs in a mixed order |
679 for (ArcIt a(_graph); a != INVALID; ++a) { |
675 int k = std::max(int(std::sqrt(double(_arc_num))), 10); |
680 _arc_id[a] = i; |
676 int i = 0, j = 0; |
681 _source[i] = _node_id[_graph.source(a)]; |
677 for (ArcIt a(_graph); a != INVALID; ++a) { |
682 _target[i] = _node_id[_graph.target(a)]; |
678 _arc_id[a] = i; |
683 if ((i += k) >= _arc_num) i = (i % k) + 1; |
679 _source[i] = _node_id[_graph.source(a)]; |
|
680 _target[i] = _node_id[_graph.target(a)]; |
|
681 if ((i += k) >= _arc_num) i = ++j; |
|
682 } |
|
683 } else { |
|
684 // Store the arcs in the original order |
|
685 int i = 0; |
|
686 for (ArcIt a(_graph); a != INVALID; ++a, ++i) { |
|
687 _arc_id[a] = i; |
|
688 _source[i] = _node_id[_graph.source(a)]; |
|
689 _target[i] = _node_id[_graph.target(a)]; |
|
690 } |
684 } |
691 } |
685 |
692 |
686 // Reset parameters |
693 // Reset parameters |
687 reset(); |
694 reset(); |
688 } |
695 } |