| ... |
... |
@@ -429,2949 +429,3009 @@
|
| 429 |
429 |
///@{
|
| 430 |
430 |
|
| 431 |
431 |
/// \brief Set the initial matching to the empty matching.
|
| 432 |
432 |
///
|
| 433 |
433 |
/// This function sets the initial matching to the empty matching.
|
| 434 |
434 |
void init() {
|
| 435 |
435 |
createStructures();
|
| 436 |
436 |
for(NodeIt n(_graph); n != INVALID; ++n) {
|
| 437 |
437 |
(*_matching)[n] = INVALID;
|
| 438 |
438 |
(*_status)[n] = UNMATCHED;
|
| 439 |
439 |
}
|
| 440 |
440 |
}
|
| 441 |
441 |
|
| 442 |
442 |
/// \brief Find an initial matching in a greedy way.
|
| 443 |
443 |
///
|
| 444 |
444 |
/// This function finds an initial matching in a greedy way.
|
| 445 |
445 |
void greedyInit() {
|
| 446 |
446 |
createStructures();
|
| 447 |
447 |
for (NodeIt n(_graph); n != INVALID; ++n) {
|
| 448 |
448 |
(*_matching)[n] = INVALID;
|
| 449 |
449 |
(*_status)[n] = UNMATCHED;
|
| 450 |
450 |
}
|
| 451 |
451 |
for (NodeIt n(_graph); n != INVALID; ++n) {
|
| 452 |
452 |
if ((*_matching)[n] == INVALID) {
|
| 453 |
453 |
for (OutArcIt a(_graph, n); a != INVALID ; ++a) {
|
| 454 |
454 |
Node v = _graph.target(a);
|
| 455 |
455 |
if ((*_matching)[v] == INVALID && v != n) {
|
| 456 |
456 |
(*_matching)[n] = a;
|
| 457 |
457 |
(*_status)[n] = MATCHED;
|
| 458 |
458 |
(*_matching)[v] = _graph.oppositeArc(a);
|
| 459 |
459 |
(*_status)[v] = MATCHED;
|
| 460 |
460 |
break;
|
| 461 |
461 |
}
|
| 462 |
462 |
}
|
| 463 |
463 |
}
|
| 464 |
464 |
}
|
| 465 |
465 |
}
|
| 466 |
466 |
|
| 467 |
467 |
|
| 468 |
468 |
/// \brief Initialize the matching from a map.
|
| 469 |
469 |
///
|
| 470 |
470 |
/// This function initializes the matching from a \c bool valued edge
|
| 471 |
471 |
/// map. This map should have the property that there are no two incident
|
| 472 |
472 |
/// edges with \c true value, i.e. it really contains a matching.
|
| 473 |
473 |
/// \return \c true if the map contains a matching.
|
| 474 |
474 |
template <typename MatchingMap>
|
| 475 |
475 |
bool matchingInit(const MatchingMap& matching) {
|
| 476 |
476 |
createStructures();
|
| 477 |
477 |
|
| 478 |
478 |
for (NodeIt n(_graph); n != INVALID; ++n) {
|
| 479 |
479 |
(*_matching)[n] = INVALID;
|
| 480 |
480 |
(*_status)[n] = UNMATCHED;
|
| 481 |
481 |
}
|
| 482 |
482 |
for(EdgeIt e(_graph); e!=INVALID; ++e) {
|
| 483 |
483 |
if (matching[e]) {
|
| 484 |
484 |
|
| 485 |
485 |
Node u = _graph.u(e);
|
| 486 |
486 |
if ((*_matching)[u] != INVALID) return false;
|
| 487 |
487 |
(*_matching)[u] = _graph.direct(e, true);
|
| 488 |
488 |
(*_status)[u] = MATCHED;
|
| 489 |
489 |
|
| 490 |
490 |
Node v = _graph.v(e);
|
| 491 |
491 |
if ((*_matching)[v] != INVALID) return false;
|
| 492 |
492 |
(*_matching)[v] = _graph.direct(e, false);
|
| 493 |
493 |
(*_status)[v] = MATCHED;
|
| 494 |
494 |
}
|
| 495 |
495 |
}
|
| 496 |
496 |
return true;
|
| 497 |
497 |
}
|
| 498 |
498 |
|
| 499 |
499 |
/// \brief Start Edmonds' algorithm
|
| 500 |
500 |
///
|
| 501 |
501 |
/// This function runs the original Edmonds' algorithm.
|
| 502 |
502 |
///
|
| 503 |
503 |
/// \pre \ref init(), \ref greedyInit() or \ref matchingInit() must be
|
| 504 |
504 |
/// called before using this function.
|
| 505 |
505 |
void startSparse() {
|
| 506 |
506 |
for(NodeIt n(_graph); n != INVALID; ++n) {
|
| 507 |
507 |
if ((*_status)[n] == UNMATCHED) {
|
| 508 |
508 |
(*_blossom_rep)[_blossom_set->insert(n)] = n;
|
| 509 |
509 |
_tree_set->insert(n);
|
| 510 |
510 |
(*_status)[n] = EVEN;
|
| 511 |
511 |
processSparse(n);
|
| 512 |
512 |
}
|
| 513 |
513 |
}
|
| 514 |
514 |
}
|
| 515 |
515 |
|
| 516 |
516 |
/// \brief Start Edmonds' algorithm with a heuristic improvement
|
| 517 |
517 |
/// for dense graphs
|
| 518 |
518 |
///
|
| 519 |
519 |
/// This function runs Edmonds' algorithm with a heuristic of postponing
|
| 520 |
520 |
/// shrinks, therefore resulting in a faster algorithm for dense graphs.
|
| 521 |
521 |
///
|
| 522 |
522 |
/// \pre \ref init(), \ref greedyInit() or \ref matchingInit() must be
|
| 523 |
523 |
/// called before using this function.
|
| 524 |
524 |
void startDense() {
|
| 525 |
525 |
for(NodeIt n(_graph); n != INVALID; ++n) {
|
| 526 |
526 |
if ((*_status)[n] == UNMATCHED) {
|
| 527 |
527 |
(*_blossom_rep)[_blossom_set->insert(n)] = n;
|
| 528 |
528 |
_tree_set->insert(n);
|
| 529 |
529 |
(*_status)[n] = EVEN;
|
| 530 |
530 |
processDense(n);
|
| 531 |
531 |
}
|
| 532 |
532 |
}
|
| 533 |
533 |
}
|
| 534 |
534 |
|
| 535 |
535 |
|
| 536 |
536 |
/// \brief Run Edmonds' algorithm
|
| 537 |
537 |
///
|
| 538 |
538 |
/// This function runs Edmonds' algorithm. An additional heuristic of
|
| 539 |
539 |
/// postponing shrinks is used for relatively dense graphs
|
| 540 |
540 |
/// (for which <tt>m>=2*n</tt> holds).
|
| 541 |
541 |
void run() {
|
| 542 |
542 |
if (countEdges(_graph) < 2 * countNodes(_graph)) {
|
| 543 |
543 |
greedyInit();
|
| 544 |
544 |
startSparse();
|
| 545 |
545 |
} else {
|
| 546 |
546 |
init();
|
| 547 |
547 |
startDense();
|
| 548 |
548 |
}
|
| 549 |
549 |
}
|
| 550 |
550 |
|
| 551 |
551 |
/// @}
|
| 552 |
552 |
|
| 553 |
553 |
/// \name Primal Solution
|
| 554 |
554 |
/// Functions to get the primal solution, i.e. the maximum matching.
|
| 555 |
555 |
|
| 556 |
556 |
/// @{
|
| 557 |
557 |
|
| 558 |
558 |
/// \brief Return the size (cardinality) of the matching.
|
| 559 |
559 |
///
|
| 560 |
560 |
/// This function returns the size (cardinality) of the current matching.
|
| 561 |
561 |
/// After run() it returns the size of the maximum matching in the graph.
|
| 562 |
562 |
int matchingSize() const {
|
| 563 |
563 |
int size = 0;
|
| 564 |
564 |
for (NodeIt n(_graph); n != INVALID; ++n) {
|
| 565 |
565 |
if ((*_matching)[n] != INVALID) {
|
| 566 |
566 |
++size;
|
| 567 |
567 |
}
|
| 568 |
568 |
}
|
| 569 |
569 |
return size / 2;
|
| 570 |
570 |
}
|
| 571 |
571 |
|
| 572 |
572 |
/// \brief Return \c true if the given edge is in the matching.
|
| 573 |
573 |
///
|
| 574 |
574 |
/// This function returns \c true if the given edge is in the current
|
| 575 |
575 |
/// matching.
|
| 576 |
576 |
bool matching(const Edge& edge) const {
|
| 577 |
577 |
return edge == (*_matching)[_graph.u(edge)];
|
| 578 |
578 |
}
|
| 579 |
579 |
|
| 580 |
580 |
/// \brief Return the matching arc (or edge) incident to the given node.
|
| 581 |
581 |
///
|
| 582 |
582 |
/// This function returns the matching arc (or edge) incident to the
|
| 583 |
583 |
/// given node in the current matching or \c INVALID if the node is
|
| 584 |
584 |
/// not covered by the matching.
|
| 585 |
585 |
Arc matching(const Node& n) const {
|
| 586 |
586 |
return (*_matching)[n];
|
| 587 |
587 |
}
|
| 588 |
588 |
|
| 589 |
589 |
/// \brief Return a const reference to the matching map.
|
| 590 |
590 |
///
|
| 591 |
591 |
/// This function returns a const reference to a node map that stores
|
| 592 |
592 |
/// the matching arc (or edge) incident to each node.
|
| 593 |
593 |
const MatchingMap& matchingMap() const {
|
| 594 |
594 |
return *_matching;
|
| 595 |
595 |
}
|
| 596 |
596 |
|
| 597 |
597 |
/// \brief Return the mate of the given node.
|
| 598 |
598 |
///
|
| 599 |
599 |
/// This function returns the mate of the given node in the current
|
| 600 |
600 |
/// matching or \c INVALID if the node is not covered by the matching.
|
| 601 |
601 |
Node mate(const Node& n) const {
|
| 602 |
602 |
return (*_matching)[n] != INVALID ?
|
| 603 |
603 |
_graph.target((*_matching)[n]) : INVALID;
|
| 604 |
604 |
}
|
| 605 |
605 |
|
| 606 |
606 |
/// @}
|
| 607 |
607 |
|
| 608 |
608 |
/// \name Dual Solution
|
| 609 |
609 |
/// Functions to get the dual solution, i.e. the Gallai-Edmonds
|
| 610 |
610 |
/// decomposition.
|
| 611 |
611 |
|
| 612 |
612 |
/// @{
|
| 613 |
613 |
|
| 614 |
614 |
/// \brief Return the status of the given node in the Edmonds-Gallai
|
| 615 |
615 |
/// decomposition.
|
| 616 |
616 |
///
|
| 617 |
617 |
/// This function returns the \ref Status "status" of the given node
|
| 618 |
618 |
/// in the Edmonds-Gallai decomposition.
|
| 619 |
619 |
Status status(const Node& n) const {
|
| 620 |
620 |
return (*_status)[n];
|
| 621 |
621 |
}
|
| 622 |
622 |
|
| 623 |
623 |
/// \brief Return a const reference to the status map, which stores
|
| 624 |
624 |
/// the Edmonds-Gallai decomposition.
|
| 625 |
625 |
///
|
| 626 |
626 |
/// This function returns a const reference to a node map that stores the
|
| 627 |
627 |
/// \ref Status "status" of each node in the Edmonds-Gallai decomposition.
|
| 628 |
628 |
const StatusMap& statusMap() const {
|
| 629 |
629 |
return *_status;
|
| 630 |
630 |
}
|
| 631 |
631 |
|
| 632 |
632 |
/// \brief Return \c true if the given node is in the barrier.
|
| 633 |
633 |
///
|
| 634 |
634 |
/// This function returns \c true if the given node is in the barrier.
|
| 635 |
635 |
bool barrier(const Node& n) const {
|
| 636 |
636 |
return (*_status)[n] == ODD;
|
| 637 |
637 |
}
|
| 638 |
638 |
|
| 639 |
639 |
/// @}
|
| 640 |
640 |
|
| 641 |
641 |
};
|
| 642 |
642 |
|
| 643 |
643 |
/// \ingroup matching
|
| 644 |
644 |
///
|
| 645 |
645 |
/// \brief Weighted matching in general graphs
|
| 646 |
646 |
///
|
| 647 |
647 |
/// This class provides an efficient implementation of Edmond's
|
| 648 |
648 |
/// maximum weighted matching algorithm. The implementation is based
|
| 649 |
649 |
/// on extensive use of priority queues and provides
|
| 650 |
650 |
/// \f$O(nm\log n)\f$ time complexity.
|
| 651 |
651 |
///
|
| 652 |
652 |
/// The maximum weighted matching problem is to find a subset of the
|
| 653 |
653 |
/// edges in an undirected graph with maximum overall weight for which
|
| 654 |
654 |
/// each node has at most one incident edge.
|
| 655 |
655 |
/// It can be formulated with the following linear program.
|
| 656 |
656 |
/// \f[ \sum_{e \in \delta(u)}x_e \le 1 \quad \forall u\in V\f]
|
| 657 |
657 |
/** \f[ \sum_{e \in \gamma(B)}x_e \le \frac{\vert B \vert - 1}{2}
|
| 658 |
658 |
\quad \forall B\in\mathcal{O}\f] */
|
| 659 |
659 |
/// \f[x_e \ge 0\quad \forall e\in E\f]
|
| 660 |
660 |
/// \f[\max \sum_{e\in E}x_ew_e\f]
|
| 661 |
661 |
/// where \f$\delta(X)\f$ is the set of edges incident to a node in
|
| 662 |
662 |
/// \f$X\f$, \f$\gamma(X)\f$ is the set of edges with both ends in
|
| 663 |
663 |
/// \f$X\f$ and \f$\mathcal{O}\f$ is the set of odd cardinality
|
| 664 |
664 |
/// subsets of the nodes.
|
| 665 |
665 |
///
|
| 666 |
666 |
/// The algorithm calculates an optimal matching and a proof of the
|
| 667 |
667 |
/// optimality. The solution of the dual problem can be used to check
|
| 668 |
668 |
/// the result of the algorithm. The dual linear problem is the
|
| 669 |
669 |
/// following.
|
| 670 |
670 |
/** \f[ y_u + y_v + \sum_{B \in \mathcal{O}, uv \in \gamma(B)}
|
| 671 |
671 |
z_B \ge w_{uv} \quad \forall uv\in E\f] */
|
| 672 |
672 |
/// \f[y_u \ge 0 \quad \forall u \in V\f]
|
| 673 |
673 |
/// \f[z_B \ge 0 \quad \forall B \in \mathcal{O}\f]
|
| 674 |
674 |
/** \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}}
|
| 675 |
675 |
\frac{\vert B \vert - 1}{2}z_B\f] */
|
| 676 |
676 |
///
|
| 677 |
677 |
/// The algorithm can be executed with the run() function.
|
| 678 |
678 |
/// After it the matching (the primal solution) and the dual solution
|
| 679 |
679 |
/// can be obtained using the query functions and the
|
| 680 |
680 |
/// \ref MaxWeightedMatching::BlossomIt "BlossomIt" nested class,
|
| 681 |
681 |
/// which is able to iterate on the nodes of a blossom.
|
| 682 |
682 |
/// If the value type is integer, then the dual solution is multiplied
|
| 683 |
683 |
/// by \ref MaxWeightedMatching::dualScale "4".
|
| 684 |
684 |
///
|
| 685 |
685 |
/// \tparam GR The undirected graph type the algorithm runs on.
|
| 686 |
686 |
/// \tparam WM The type edge weight map. The default type is
|
| 687 |
687 |
/// \ref concepts::Graph::EdgeMap "GR::EdgeMap<int>".
|
| 688 |
688 |
#ifdef DOXYGEN
|
| 689 |
689 |
template <typename GR, typename WM>
|
| 690 |
690 |
#else
|
| 691 |
691 |
template <typename GR,
|
| 692 |
692 |
typename WM = typename GR::template EdgeMap<int> >
|
| 693 |
693 |
#endif
|
| 694 |
694 |
class MaxWeightedMatching {
|
| 695 |
695 |
public:
|
| 696 |
696 |
|
| 697 |
697 |
/// The graph type of the algorithm
|
| 698 |
698 |
typedef GR Graph;
|
| 699 |
699 |
/// The type of the edge weight map
|
| 700 |
700 |
typedef WM WeightMap;
|
| 701 |
701 |
/// The value type of the edge weights
|
| 702 |
702 |
typedef typename WeightMap::Value Value;
|
| 703 |
703 |
|
| 704 |
704 |
/// The type of the matching map
|
| 705 |
705 |
typedef typename Graph::template NodeMap<typename Graph::Arc>
|
| 706 |
706 |
MatchingMap;
|
| 707 |
707 |
|
| 708 |
708 |
/// \brief Scaling factor for dual solution
|
| 709 |
709 |
///
|
| 710 |
710 |
/// Scaling factor for dual solution. It is equal to 4 or 1
|
| 711 |
711 |
/// according to the value type.
|
| 712 |
712 |
static const int dualScale =
|
| 713 |
713 |
std::numeric_limits<Value>::is_integer ? 4 : 1;
|
| 714 |
714 |
|
| 715 |
715 |
private:
|
| 716 |
716 |
|
| 717 |
717 |
TEMPLATE_GRAPH_TYPEDEFS(Graph);
|
| 718 |
718 |
|
| 719 |
719 |
typedef typename Graph::template NodeMap<Value> NodePotential;
|
| 720 |
720 |
typedef std::vector<Node> BlossomNodeList;
|
| 721 |
721 |
|
| 722 |
722 |
struct BlossomVariable {
|
| 723 |
723 |
int begin, end;
|
| 724 |
724 |
Value value;
|
| 725 |
725 |
|
| 726 |
726 |
BlossomVariable(int _begin, int _end, Value _value)
|
| 727 |
727 |
: begin(_begin), end(_end), value(_value) {}
|
| 728 |
728 |
|
| 729 |
729 |
};
|
| 730 |
730 |
|
| 731 |
731 |
typedef std::vector<BlossomVariable> BlossomPotential;
|
| 732 |
732 |
|
| 733 |
733 |
const Graph& _graph;
|
| 734 |
734 |
const WeightMap& _weight;
|
| 735 |
735 |
|
| 736 |
736 |
MatchingMap* _matching;
|
| 737 |
737 |
|
| 738 |
738 |
NodePotential* _node_potential;
|
| 739 |
739 |
|
| 740 |
740 |
BlossomPotential _blossom_potential;
|
| 741 |
741 |
BlossomNodeList _blossom_node_list;
|
| 742 |
742 |
|
| 743 |
743 |
int _node_num;
|
| 744 |
744 |
int _blossom_num;
|
| 745 |
745 |
|
| 746 |
746 |
typedef RangeMap<int> IntIntMap;
|
| 747 |
747 |
|
| 748 |
748 |
enum Status {
|
| 749 |
749 |
EVEN = -1, MATCHED = 0, ODD = 1
|
| 750 |
750 |
};
|
| 751 |
751 |
|
| 752 |
752 |
typedef HeapUnionFind<Value, IntNodeMap> BlossomSet;
|
| 753 |
753 |
struct BlossomData {
|
| 754 |
754 |
int tree;
|
| 755 |
755 |
Status status;
|
| 756 |
756 |
Arc pred, next;
|
| 757 |
757 |
Value pot, offset;
|
| 758 |
758 |
Node base;
|
| 759 |
759 |
};
|
| 760 |
760 |
|
| 761 |
761 |
IntNodeMap *_blossom_index;
|
| 762 |
762 |
BlossomSet *_blossom_set;
|
| 763 |
763 |
RangeMap<BlossomData>* _blossom_data;
|
| 764 |
764 |
|
| 765 |
765 |
IntNodeMap *_node_index;
|
| 766 |
766 |
IntArcMap *_node_heap_index;
|
| 767 |
767 |
|
| 768 |
768 |
struct NodeData {
|
| 769 |
769 |
|
| 770 |
770 |
NodeData(IntArcMap& node_heap_index)
|
| 771 |
771 |
: heap(node_heap_index) {}
|
| 772 |
772 |
|
| 773 |
773 |
int blossom;
|
| 774 |
774 |
Value pot;
|
| 775 |
775 |
BinHeap<Value, IntArcMap> heap;
|
| 776 |
776 |
std::map<int, Arc> heap_index;
|
| 777 |
777 |
|
| 778 |
778 |
int tree;
|
| 779 |
779 |
};
|
| 780 |
780 |
|
| 781 |
781 |
RangeMap<NodeData>* _node_data;
|
| 782 |
782 |
|
| 783 |
783 |
typedef ExtendFindEnum<IntIntMap> TreeSet;
|
| 784 |
784 |
|
| 785 |
785 |
IntIntMap *_tree_set_index;
|
| 786 |
786 |
TreeSet *_tree_set;
|
| 787 |
787 |
|
| 788 |
788 |
IntNodeMap *_delta1_index;
|
| 789 |
789 |
BinHeap<Value, IntNodeMap> *_delta1;
|
| 790 |
790 |
|
| 791 |
791 |
IntIntMap *_delta2_index;
|
| 792 |
792 |
BinHeap<Value, IntIntMap> *_delta2;
|
| 793 |
793 |
|
| 794 |
794 |
IntEdgeMap *_delta3_index;
|
| 795 |
795 |
BinHeap<Value, IntEdgeMap> *_delta3;
|
| 796 |
796 |
|
| 797 |
797 |
IntIntMap *_delta4_index;
|
| 798 |
798 |
BinHeap<Value, IntIntMap> *_delta4;
|
| 799 |
799 |
|
| 800 |
800 |
Value _delta_sum;
|
| 801 |
801 |
int _unmatched;
|
| 802 |
802 |
|
| 803 |
803 |
typedef MaxWeightedFractionalMatching<Graph, WeightMap> FractionalMatching;
|
| 804 |
804 |
FractionalMatching *_fractional;
|
| 805 |
805 |
|
| 806 |
806 |
void createStructures() {
|
| 807 |
807 |
_node_num = countNodes(_graph);
|
| 808 |
808 |
_blossom_num = _node_num * 3 / 2;
|
| 809 |
809 |
|
| 810 |
810 |
if (!_matching) {
|
| 811 |
811 |
_matching = new MatchingMap(_graph);
|
| 812 |
812 |
}
|
|
813 |
|
| 813 |
814 |
if (!_node_potential) {
|
| 814 |
815 |
_node_potential = new NodePotential(_graph);
|
| 815 |
816 |
}
|
|
817 |
|
| 816 |
818 |
if (!_blossom_set) {
|
| 817 |
819 |
_blossom_index = new IntNodeMap(_graph);
|
| 818 |
820 |
_blossom_set = new BlossomSet(*_blossom_index);
|
| 819 |
821 |
_blossom_data = new RangeMap<BlossomData>(_blossom_num);
|
|
822 |
} else if (_blossom_data->size() != _blossom_num) {
|
|
823 |
delete _blossom_data;
|
|
824 |
_blossom_data = new RangeMap<BlossomData>(_blossom_num);
|
| 820 |
825 |
}
|
| 821 |
826 |
|
| 822 |
827 |
if (!_node_index) {
|
| 823 |
828 |
_node_index = new IntNodeMap(_graph);
|
| 824 |
829 |
_node_heap_index = new IntArcMap(_graph);
|
| 825 |
830 |
_node_data = new RangeMap<NodeData>(_node_num,
|
| 826 |
|
NodeData(*_node_heap_index));
|
|
831 |
NodeData(*_node_heap_index));
|
|
832 |
} else {
|
|
833 |
delete _node_data;
|
|
834 |
_node_data = new RangeMap<NodeData>(_node_num,
|
|
835 |
NodeData(*_node_heap_index));
|
| 827 |
836 |
}
|
| 828 |
837 |
|
| 829 |
838 |
if (!_tree_set) {
|
| 830 |
839 |
_tree_set_index = new IntIntMap(_blossom_num);
|
| 831 |
840 |
_tree_set = new TreeSet(*_tree_set_index);
|
|
841 |
} else {
|
|
842 |
_tree_set_index->resize(_blossom_num);
|
| 832 |
843 |
}
|
|
844 |
|
| 833 |
845 |
if (!_delta1) {
|
| 834 |
846 |
_delta1_index = new IntNodeMap(_graph);
|
| 835 |
847 |
_delta1 = new BinHeap<Value, IntNodeMap>(*_delta1_index);
|
| 836 |
848 |
}
|
|
849 |
|
| 837 |
850 |
if (!_delta2) {
|
| 838 |
851 |
_delta2_index = new IntIntMap(_blossom_num);
|
| 839 |
852 |
_delta2 = new BinHeap<Value, IntIntMap>(*_delta2_index);
|
|
853 |
} else {
|
|
854 |
_delta2_index->resize(_blossom_num);
|
| 840 |
855 |
}
|
|
856 |
|
| 841 |
857 |
if (!_delta3) {
|
| 842 |
858 |
_delta3_index = new IntEdgeMap(_graph);
|
| 843 |
859 |
_delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index);
|
| 844 |
860 |
}
|
|
861 |
|
| 845 |
862 |
if (!_delta4) {
|
| 846 |
863 |
_delta4_index = new IntIntMap(_blossom_num);
|
| 847 |
864 |
_delta4 = new BinHeap<Value, IntIntMap>(*_delta4_index);
|
|
865 |
} else {
|
|
866 |
_delta4_index->resize(_blossom_num);
|
| 848 |
867 |
}
|
| 849 |
868 |
}
|
| 850 |
869 |
|
| 851 |
870 |
void destroyStructures() {
|
| 852 |
871 |
if (_matching) {
|
| 853 |
872 |
delete _matching;
|
| 854 |
873 |
}
|
| 855 |
874 |
if (_node_potential) {
|
| 856 |
875 |
delete _node_potential;
|
| 857 |
876 |
}
|
| 858 |
877 |
if (_blossom_set) {
|
| 859 |
878 |
delete _blossom_index;
|
| 860 |
879 |
delete _blossom_set;
|
| 861 |
880 |
delete _blossom_data;
|
| 862 |
881 |
}
|
| 863 |
882 |
|
| 864 |
883 |
if (_node_index) {
|
| 865 |
884 |
delete _node_index;
|
| 866 |
885 |
delete _node_heap_index;
|
| 867 |
886 |
delete _node_data;
|
| 868 |
887 |
}
|
| 869 |
888 |
|
| 870 |
889 |
if (_tree_set) {
|
| 871 |
890 |
delete _tree_set_index;
|
| 872 |
891 |
delete _tree_set;
|
| 873 |
892 |
}
|
| 874 |
893 |
if (_delta1) {
|
| 875 |
894 |
delete _delta1_index;
|
| 876 |
895 |
delete _delta1;
|
| 877 |
896 |
}
|
| 878 |
897 |
if (_delta2) {
|
| 879 |
898 |
delete _delta2_index;
|
| 880 |
899 |
delete _delta2;
|
| 881 |
900 |
}
|
| 882 |
901 |
if (_delta3) {
|
| 883 |
902 |
delete _delta3_index;
|
| 884 |
903 |
delete _delta3;
|
| 885 |
904 |
}
|
| 886 |
905 |
if (_delta4) {
|
| 887 |
906 |
delete _delta4_index;
|
| 888 |
907 |
delete _delta4;
|
| 889 |
908 |
}
|
| 890 |
909 |
}
|
| 891 |
910 |
|
| 892 |
911 |
void matchedToEven(int blossom, int tree) {
|
| 893 |
912 |
if (_delta2->state(blossom) == _delta2->IN_HEAP) {
|
| 894 |
913 |
_delta2->erase(blossom);
|
| 895 |
914 |
}
|
| 896 |
915 |
|
| 897 |
916 |
if (!_blossom_set->trivial(blossom)) {
|
| 898 |
917 |
(*_blossom_data)[blossom].pot -=
|
| 899 |
918 |
2 * (_delta_sum - (*_blossom_data)[blossom].offset);
|
| 900 |
919 |
}
|
| 901 |
920 |
|
| 902 |
921 |
for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
|
| 903 |
922 |
n != INVALID; ++n) {
|
| 904 |
923 |
|
| 905 |
924 |
_blossom_set->increase(n, std::numeric_limits<Value>::max());
|
| 906 |
925 |
int ni = (*_node_index)[n];
|
| 907 |
926 |
|
| 908 |
927 |
(*_node_data)[ni].heap.clear();
|
| 909 |
928 |
(*_node_data)[ni].heap_index.clear();
|
| 910 |
929 |
|
| 911 |
930 |
(*_node_data)[ni].pot += _delta_sum - (*_blossom_data)[blossom].offset;
|
| 912 |
931 |
|
| 913 |
932 |
_delta1->push(n, (*_node_data)[ni].pot);
|
| 914 |
933 |
|
| 915 |
934 |
for (InArcIt e(_graph, n); e != INVALID; ++e) {
|
| 916 |
935 |
Node v = _graph.source(e);
|
| 917 |
936 |
int vb = _blossom_set->find(v);
|
| 918 |
937 |
int vi = (*_node_index)[v];
|
| 919 |
938 |
|
| 920 |
939 |
Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
|
| 921 |
940 |
dualScale * _weight[e];
|
| 922 |
941 |
|
| 923 |
942 |
if ((*_blossom_data)[vb].status == EVEN) {
|
| 924 |
943 |
if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
|
| 925 |
944 |
_delta3->push(e, rw / 2);
|
| 926 |
945 |
}
|
| 927 |
946 |
} else {
|
| 928 |
947 |
typename std::map<int, Arc>::iterator it =
|
| 929 |
948 |
(*_node_data)[vi].heap_index.find(tree);
|
| 930 |
949 |
|
| 931 |
950 |
if (it != (*_node_data)[vi].heap_index.end()) {
|
| 932 |
951 |
if ((*_node_data)[vi].heap[it->second] > rw) {
|
| 933 |
952 |
(*_node_data)[vi].heap.replace(it->second, e);
|
| 934 |
953 |
(*_node_data)[vi].heap.decrease(e, rw);
|
| 935 |
954 |
it->second = e;
|
| 936 |
955 |
}
|
| 937 |
956 |
} else {
|
| 938 |
957 |
(*_node_data)[vi].heap.push(e, rw);
|
| 939 |
958 |
(*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
|
| 940 |
959 |
}
|
| 941 |
960 |
|
| 942 |
961 |
if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
|
| 943 |
962 |
_blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
|
| 944 |
963 |
|
| 945 |
964 |
if ((*_blossom_data)[vb].status == MATCHED) {
|
| 946 |
965 |
if (_delta2->state(vb) != _delta2->IN_HEAP) {
|
| 947 |
966 |
_delta2->push(vb, _blossom_set->classPrio(vb) -
|
| 948 |
967 |
(*_blossom_data)[vb].offset);
|
| 949 |
968 |
} else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
|
| 950 |
969 |
(*_blossom_data)[vb].offset) {
|
| 951 |
970 |
_delta2->decrease(vb, _blossom_set->classPrio(vb) -
|
| 952 |
971 |
(*_blossom_data)[vb].offset);
|
| 953 |
972 |
}
|
| 954 |
973 |
}
|
| 955 |
974 |
}
|
| 956 |
975 |
}
|
| 957 |
976 |
}
|
| 958 |
977 |
}
|
| 959 |
978 |
(*_blossom_data)[blossom].offset = 0;
|
| 960 |
979 |
}
|
| 961 |
980 |
|
| 962 |
981 |
void matchedToOdd(int blossom) {
|
| 963 |
982 |
if (_delta2->state(blossom) == _delta2->IN_HEAP) {
|
| 964 |
983 |
_delta2->erase(blossom);
|
| 965 |
984 |
}
|
| 966 |
985 |
(*_blossom_data)[blossom].offset += _delta_sum;
|
| 967 |
986 |
if (!_blossom_set->trivial(blossom)) {
|
| 968 |
987 |
_delta4->push(blossom, (*_blossom_data)[blossom].pot / 2 +
|
| 969 |
988 |
(*_blossom_data)[blossom].offset);
|
| 970 |
989 |
}
|
| 971 |
990 |
}
|
| 972 |
991 |
|
| 973 |
992 |
void evenToMatched(int blossom, int tree) {
|
| 974 |
993 |
if (!_blossom_set->trivial(blossom)) {
|
| 975 |
994 |
(*_blossom_data)[blossom].pot += 2 * _delta_sum;
|
| 976 |
995 |
}
|
| 977 |
996 |
|
| 978 |
997 |
for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
|
| 979 |
998 |
n != INVALID; ++n) {
|
| 980 |
999 |
int ni = (*_node_index)[n];
|
| 981 |
1000 |
(*_node_data)[ni].pot -= _delta_sum;
|
| 982 |
1001 |
|
| 983 |
1002 |
_delta1->erase(n);
|
| 984 |
1003 |
|
| 985 |
1004 |
for (InArcIt e(_graph, n); e != INVALID; ++e) {
|
| 986 |
1005 |
Node v = _graph.source(e);
|
| 987 |
1006 |
int vb = _blossom_set->find(v);
|
| 988 |
1007 |
int vi = (*_node_index)[v];
|
| 989 |
1008 |
|
| 990 |
1009 |
Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
|
| 991 |
1010 |
dualScale * _weight[e];
|
| 992 |
1011 |
|
| 993 |
1012 |
if (vb == blossom) {
|
| 994 |
1013 |
if (_delta3->state(e) == _delta3->IN_HEAP) {
|
| 995 |
1014 |
_delta3->erase(e);
|
| 996 |
1015 |
}
|
| 997 |
1016 |
} else if ((*_blossom_data)[vb].status == EVEN) {
|
| 998 |
1017 |
|
| 999 |
1018 |
if (_delta3->state(e) == _delta3->IN_HEAP) {
|
| 1000 |
1019 |
_delta3->erase(e);
|
| 1001 |
1020 |
}
|
| 1002 |
1021 |
|
| 1003 |
1022 |
int vt = _tree_set->find(vb);
|
| 1004 |
1023 |
|
| 1005 |
1024 |
if (vt != tree) {
|
| 1006 |
1025 |
|
| 1007 |
1026 |
Arc r = _graph.oppositeArc(e);
|
| 1008 |
1027 |
|
| 1009 |
1028 |
typename std::map<int, Arc>::iterator it =
|
| 1010 |
1029 |
(*_node_data)[ni].heap_index.find(vt);
|
| 1011 |
1030 |
|
| 1012 |
1031 |
if (it != (*_node_data)[ni].heap_index.end()) {
|
| 1013 |
1032 |
if ((*_node_data)[ni].heap[it->second] > rw) {
|
| 1014 |
1033 |
(*_node_data)[ni].heap.replace(it->second, r);
|
| 1015 |
1034 |
(*_node_data)[ni].heap.decrease(r, rw);
|
| 1016 |
1035 |
it->second = r;
|
| 1017 |
1036 |
}
|
| 1018 |
1037 |
} else {
|
| 1019 |
1038 |
(*_node_data)[ni].heap.push(r, rw);
|
| 1020 |
1039 |
(*_node_data)[ni].heap_index.insert(std::make_pair(vt, r));
|
| 1021 |
1040 |
}
|
| 1022 |
1041 |
|
| 1023 |
1042 |
if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) {
|
| 1024 |
1043 |
_blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
|
| 1025 |
1044 |
|
| 1026 |
1045 |
if (_delta2->state(blossom) != _delta2->IN_HEAP) {
|
| 1027 |
1046 |
_delta2->push(blossom, _blossom_set->classPrio(blossom) -
|
| 1028 |
1047 |
(*_blossom_data)[blossom].offset);
|
| 1029 |
1048 |
} else if ((*_delta2)[blossom] >
|
| 1030 |
1049 |
_blossom_set->classPrio(blossom) -
|
| 1031 |
1050 |
(*_blossom_data)[blossom].offset){
|
| 1032 |
1051 |
_delta2->decrease(blossom, _blossom_set->classPrio(blossom) -
|
| 1033 |
1052 |
(*_blossom_data)[blossom].offset);
|
| 1034 |
1053 |
}
|
| 1035 |
1054 |
}
|
| 1036 |
1055 |
}
|
| 1037 |
1056 |
} else {
|
| 1038 |
1057 |
|
| 1039 |
1058 |
typename std::map<int, Arc>::iterator it =
|
| 1040 |
1059 |
(*_node_data)[vi].heap_index.find(tree);
|
| 1041 |
1060 |
|
| 1042 |
1061 |
if (it != (*_node_data)[vi].heap_index.end()) {
|
| 1043 |
1062 |
(*_node_data)[vi].heap.erase(it->second);
|
| 1044 |
1063 |
(*_node_data)[vi].heap_index.erase(it);
|
| 1045 |
1064 |
if ((*_node_data)[vi].heap.empty()) {
|
| 1046 |
1065 |
_blossom_set->increase(v, std::numeric_limits<Value>::max());
|
| 1047 |
1066 |
} else if ((*_blossom_set)[v] < (*_node_data)[vi].heap.prio()) {
|
| 1048 |
1067 |
_blossom_set->increase(v, (*_node_data)[vi].heap.prio());
|
| 1049 |
1068 |
}
|
| 1050 |
1069 |
|
| 1051 |
1070 |
if ((*_blossom_data)[vb].status == MATCHED) {
|
| 1052 |
1071 |
if (_blossom_set->classPrio(vb) ==
|
| 1053 |
1072 |
std::numeric_limits<Value>::max()) {
|
| 1054 |
1073 |
_delta2->erase(vb);
|
| 1055 |
1074 |
} else if ((*_delta2)[vb] < _blossom_set->classPrio(vb) -
|
| 1056 |
1075 |
(*_blossom_data)[vb].offset) {
|
| 1057 |
1076 |
_delta2->increase(vb, _blossom_set->classPrio(vb) -
|
| 1058 |
1077 |
(*_blossom_data)[vb].offset);
|
| 1059 |
1078 |
}
|
| 1060 |
1079 |
}
|
| 1061 |
1080 |
}
|
| 1062 |
1081 |
}
|
| 1063 |
1082 |
}
|
| 1064 |
1083 |
}
|
| 1065 |
1084 |
}
|
| 1066 |
1085 |
|
| 1067 |
1086 |
void oddToMatched(int blossom) {
|
| 1068 |
1087 |
(*_blossom_data)[blossom].offset -= _delta_sum;
|
| 1069 |
1088 |
|
| 1070 |
1089 |
if (_blossom_set->classPrio(blossom) !=
|
| 1071 |
1090 |
std::numeric_limits<Value>::max()) {
|
| 1072 |
1091 |
_delta2->push(blossom, _blossom_set->classPrio(blossom) -
|
| 1073 |
1092 |
(*_blossom_data)[blossom].offset);
|
| 1074 |
1093 |
}
|
| 1075 |
1094 |
|
| 1076 |
1095 |
if (!_blossom_set->trivial(blossom)) {
|
| 1077 |
1096 |
_delta4->erase(blossom);
|
| 1078 |
1097 |
}
|
| 1079 |
1098 |
}
|
| 1080 |
1099 |
|
| 1081 |
1100 |
void oddToEven(int blossom, int tree) {
|
| 1082 |
1101 |
if (!_blossom_set->trivial(blossom)) {
|
| 1083 |
1102 |
_delta4->erase(blossom);
|
| 1084 |
1103 |
(*_blossom_data)[blossom].pot -=
|
| 1085 |
1104 |
2 * (2 * _delta_sum - (*_blossom_data)[blossom].offset);
|
| 1086 |
1105 |
}
|
| 1087 |
1106 |
|
| 1088 |
1107 |
for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
|
| 1089 |
1108 |
n != INVALID; ++n) {
|
| 1090 |
1109 |
int ni = (*_node_index)[n];
|
| 1091 |
1110 |
|
| 1092 |
1111 |
_blossom_set->increase(n, std::numeric_limits<Value>::max());
|
| 1093 |
1112 |
|
| 1094 |
1113 |
(*_node_data)[ni].heap.clear();
|
| 1095 |
1114 |
(*_node_data)[ni].heap_index.clear();
|
| 1096 |
1115 |
(*_node_data)[ni].pot +=
|
| 1097 |
1116 |
2 * _delta_sum - (*_blossom_data)[blossom].offset;
|
| 1098 |
1117 |
|
| 1099 |
1118 |
_delta1->push(n, (*_node_data)[ni].pot);
|
| 1100 |
1119 |
|
| 1101 |
1120 |
for (InArcIt e(_graph, n); e != INVALID; ++e) {
|
| 1102 |
1121 |
Node v = _graph.source(e);
|
| 1103 |
1122 |
int vb = _blossom_set->find(v);
|
| 1104 |
1123 |
int vi = (*_node_index)[v];
|
| 1105 |
1124 |
|
| 1106 |
1125 |
Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
|
| 1107 |
1126 |
dualScale * _weight[e];
|
| 1108 |
1127 |
|
| 1109 |
1128 |
if ((*_blossom_data)[vb].status == EVEN) {
|
| 1110 |
1129 |
if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
|
| 1111 |
1130 |
_delta3->push(e, rw / 2);
|
| 1112 |
1131 |
}
|
| 1113 |
1132 |
} else {
|
| 1114 |
1133 |
|
| 1115 |
1134 |
typename std::map<int, Arc>::iterator it =
|
| 1116 |
1135 |
(*_node_data)[vi].heap_index.find(tree);
|
| 1117 |
1136 |
|
| 1118 |
1137 |
if (it != (*_node_data)[vi].heap_index.end()) {
|
| 1119 |
1138 |
if ((*_node_data)[vi].heap[it->second] > rw) {
|
| 1120 |
1139 |
(*_node_data)[vi].heap.replace(it->second, e);
|
| 1121 |
1140 |
(*_node_data)[vi].heap.decrease(e, rw);
|
| 1122 |
1141 |
it->second = e;
|
| 1123 |
1142 |
}
|
| 1124 |
1143 |
} else {
|
| 1125 |
1144 |
(*_node_data)[vi].heap.push(e, rw);
|
| 1126 |
1145 |
(*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
|
| 1127 |
1146 |
}
|
| 1128 |
1147 |
|
| 1129 |
1148 |
if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
|
| 1130 |
1149 |
_blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
|
| 1131 |
1150 |
|
| 1132 |
1151 |
if ((*_blossom_data)[vb].status == MATCHED) {
|
| 1133 |
1152 |
if (_delta2->state(vb) != _delta2->IN_HEAP) {
|
| 1134 |
1153 |
_delta2->push(vb, _blossom_set->classPrio(vb) -
|
| 1135 |
1154 |
(*_blossom_data)[vb].offset);
|
| 1136 |
1155 |
} else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
|
| 1137 |
1156 |
(*_blossom_data)[vb].offset) {
|
| 1138 |
1157 |
_delta2->decrease(vb, _blossom_set->classPrio(vb) -
|
| 1139 |
1158 |
(*_blossom_data)[vb].offset);
|
| 1140 |
1159 |
}
|
| 1141 |
1160 |
}
|
| 1142 |
1161 |
}
|
| 1143 |
1162 |
}
|
| 1144 |
1163 |
}
|
| 1145 |
1164 |
}
|
| 1146 |
1165 |
(*_blossom_data)[blossom].offset = 0;
|
| 1147 |
1166 |
}
|
| 1148 |
1167 |
|
| 1149 |
1168 |
void alternatePath(int even, int tree) {
|
| 1150 |
1169 |
int odd;
|
| 1151 |
1170 |
|
| 1152 |
1171 |
evenToMatched(even, tree);
|
| 1153 |
1172 |
(*_blossom_data)[even].status = MATCHED;
|
| 1154 |
1173 |
|
| 1155 |
1174 |
while ((*_blossom_data)[even].pred != INVALID) {
|
| 1156 |
1175 |
odd = _blossom_set->find(_graph.target((*_blossom_data)[even].pred));
|
| 1157 |
1176 |
(*_blossom_data)[odd].status = MATCHED;
|
| 1158 |
1177 |
oddToMatched(odd);
|
| 1159 |
1178 |
(*_blossom_data)[odd].next = (*_blossom_data)[odd].pred;
|
| 1160 |
1179 |
|
| 1161 |
1180 |
even = _blossom_set->find(_graph.target((*_blossom_data)[odd].pred));
|
| 1162 |
1181 |
(*_blossom_data)[even].status = MATCHED;
|
| 1163 |
1182 |
evenToMatched(even, tree);
|
| 1164 |
1183 |
(*_blossom_data)[even].next =
|
| 1165 |
1184 |
_graph.oppositeArc((*_blossom_data)[odd].pred);
|
| 1166 |
1185 |
}
|
| 1167 |
1186 |
|
| 1168 |
1187 |
}
|
| 1169 |
1188 |
|
| 1170 |
1189 |
void destroyTree(int tree) {
|
| 1171 |
1190 |
for (TreeSet::ItemIt b(*_tree_set, tree); b != INVALID; ++b) {
|
| 1172 |
1191 |
if ((*_blossom_data)[b].status == EVEN) {
|
| 1173 |
1192 |
(*_blossom_data)[b].status = MATCHED;
|
| 1174 |
1193 |
evenToMatched(b, tree);
|
| 1175 |
1194 |
} else if ((*_blossom_data)[b].status == ODD) {
|
| 1176 |
1195 |
(*_blossom_data)[b].status = MATCHED;
|
| 1177 |
1196 |
oddToMatched(b);
|
| 1178 |
1197 |
}
|
| 1179 |
1198 |
}
|
| 1180 |
1199 |
_tree_set->eraseClass(tree);
|
| 1181 |
1200 |
}
|
| 1182 |
1201 |
|
| 1183 |
1202 |
|
| 1184 |
1203 |
void unmatchNode(const Node& node) {
|
| 1185 |
1204 |
int blossom = _blossom_set->find(node);
|
| 1186 |
1205 |
int tree = _tree_set->find(blossom);
|
| 1187 |
1206 |
|
| 1188 |
1207 |
alternatePath(blossom, tree);
|
| 1189 |
1208 |
destroyTree(tree);
|
| 1190 |
1209 |
|
| 1191 |
1210 |
(*_blossom_data)[blossom].base = node;
|
| 1192 |
1211 |
(*_blossom_data)[blossom].next = INVALID;
|
| 1193 |
1212 |
}
|
| 1194 |
1213 |
|
| 1195 |
1214 |
void augmentOnEdge(const Edge& edge) {
|
| 1196 |
1215 |
|
| 1197 |
1216 |
int left = _blossom_set->find(_graph.u(edge));
|
| 1198 |
1217 |
int right = _blossom_set->find(_graph.v(edge));
|
| 1199 |
1218 |
|
| 1200 |
1219 |
int left_tree = _tree_set->find(left);
|
| 1201 |
1220 |
alternatePath(left, left_tree);
|
| 1202 |
1221 |
destroyTree(left_tree);
|
| 1203 |
1222 |
|
| 1204 |
1223 |
int right_tree = _tree_set->find(right);
|
| 1205 |
1224 |
alternatePath(right, right_tree);
|
| 1206 |
1225 |
destroyTree(right_tree);
|
| 1207 |
1226 |
|
| 1208 |
1227 |
(*_blossom_data)[left].next = _graph.direct(edge, true);
|
| 1209 |
1228 |
(*_blossom_data)[right].next = _graph.direct(edge, false);
|
| 1210 |
1229 |
}
|
| 1211 |
1230 |
|
| 1212 |
1231 |
void augmentOnArc(const Arc& arc) {
|
| 1213 |
1232 |
|
| 1214 |
1233 |
int left = _blossom_set->find(_graph.source(arc));
|
| 1215 |
1234 |
int right = _blossom_set->find(_graph.target(arc));
|
| 1216 |
1235 |
|
| 1217 |
1236 |
(*_blossom_data)[left].status = MATCHED;
|
| 1218 |
1237 |
|
| 1219 |
1238 |
int right_tree = _tree_set->find(right);
|
| 1220 |
1239 |
alternatePath(right, right_tree);
|
| 1221 |
1240 |
destroyTree(right_tree);
|
| 1222 |
1241 |
|
| 1223 |
1242 |
(*_blossom_data)[left].next = arc;
|
| 1224 |
1243 |
(*_blossom_data)[right].next = _graph.oppositeArc(arc);
|
| 1225 |
1244 |
}
|
| 1226 |
1245 |
|
| 1227 |
1246 |
void extendOnArc(const Arc& arc) {
|
| 1228 |
1247 |
int base = _blossom_set->find(_graph.target(arc));
|
| 1229 |
1248 |
int tree = _tree_set->find(base);
|
| 1230 |
1249 |
|
| 1231 |
1250 |
int odd = _blossom_set->find(_graph.source(arc));
|
| 1232 |
1251 |
_tree_set->insert(odd, tree);
|
| 1233 |
1252 |
(*_blossom_data)[odd].status = ODD;
|
| 1234 |
1253 |
matchedToOdd(odd);
|
| 1235 |
1254 |
(*_blossom_data)[odd].pred = arc;
|
| 1236 |
1255 |
|
| 1237 |
1256 |
int even = _blossom_set->find(_graph.target((*_blossom_data)[odd].next));
|
| 1238 |
1257 |
(*_blossom_data)[even].pred = (*_blossom_data)[even].next;
|
| 1239 |
1258 |
_tree_set->insert(even, tree);
|
| 1240 |
1259 |
(*_blossom_data)[even].status = EVEN;
|
| 1241 |
1260 |
matchedToEven(even, tree);
|
| 1242 |
1261 |
}
|
| 1243 |
1262 |
|
| 1244 |
1263 |
void shrinkOnEdge(const Edge& edge, int tree) {
|
| 1245 |
1264 |
int nca = -1;
|
| 1246 |
1265 |
std::vector<int> left_path, right_path;
|
| 1247 |
1266 |
|
| 1248 |
1267 |
{
|
| 1249 |
1268 |
std::set<int> left_set, right_set;
|
| 1250 |
1269 |
int left = _blossom_set->find(_graph.u(edge));
|
| 1251 |
1270 |
left_path.push_back(left);
|
| 1252 |
1271 |
left_set.insert(left);
|
| 1253 |
1272 |
|
| 1254 |
1273 |
int right = _blossom_set->find(_graph.v(edge));
|
| 1255 |
1274 |
right_path.push_back(right);
|
| 1256 |
1275 |
right_set.insert(right);
|
| 1257 |
1276 |
|
| 1258 |
1277 |
while (true) {
|
| 1259 |
1278 |
|
| 1260 |
1279 |
if ((*_blossom_data)[left].pred == INVALID) break;
|
| 1261 |
1280 |
|
| 1262 |
1281 |
left =
|
| 1263 |
1282 |
_blossom_set->find(_graph.target((*_blossom_data)[left].pred));
|
| 1264 |
1283 |
left_path.push_back(left);
|
| 1265 |
1284 |
left =
|
| 1266 |
1285 |
_blossom_set->find(_graph.target((*_blossom_data)[left].pred));
|
| 1267 |
1286 |
left_path.push_back(left);
|
| 1268 |
1287 |
|
| 1269 |
1288 |
left_set.insert(left);
|
| 1270 |
1289 |
|
| 1271 |
1290 |
if (right_set.find(left) != right_set.end()) {
|
| 1272 |
1291 |
nca = left;
|
| 1273 |
1292 |
break;
|
| 1274 |
1293 |
}
|
| 1275 |
1294 |
|
| 1276 |
1295 |
if ((*_blossom_data)[right].pred == INVALID) break;
|
| 1277 |
1296 |
|
| 1278 |
1297 |
right =
|
| 1279 |
1298 |
_blossom_set->find(_graph.target((*_blossom_data)[right].pred));
|
| 1280 |
1299 |
right_path.push_back(right);
|
| 1281 |
1300 |
right =
|
| 1282 |
1301 |
_blossom_set->find(_graph.target((*_blossom_data)[right].pred));
|
| 1283 |
1302 |
right_path.push_back(right);
|
| 1284 |
1303 |
|
| 1285 |
1304 |
right_set.insert(right);
|
| 1286 |
1305 |
|
| 1287 |
1306 |
if (left_set.find(right) != left_set.end()) {
|
| 1288 |
1307 |
nca = right;
|
| 1289 |
1308 |
break;
|
| 1290 |
1309 |
}
|
| 1291 |
1310 |
|
| 1292 |
1311 |
}
|
| 1293 |
1312 |
|
| 1294 |
1313 |
if (nca == -1) {
|
| 1295 |
1314 |
if ((*_blossom_data)[left].pred == INVALID) {
|
| 1296 |
1315 |
nca = right;
|
| 1297 |
1316 |
while (left_set.find(nca) == left_set.end()) {
|
| 1298 |
1317 |
nca =
|
| 1299 |
1318 |
_blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
|
| 1300 |
1319 |
right_path.push_back(nca);
|
| 1301 |
1320 |
nca =
|
| 1302 |
1321 |
_blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
|
| 1303 |
1322 |
right_path.push_back(nca);
|
| 1304 |
1323 |
}
|
| 1305 |
1324 |
} else {
|
| 1306 |
1325 |
nca = left;
|
| 1307 |
1326 |
while (right_set.find(nca) == right_set.end()) {
|
| 1308 |
1327 |
nca =
|
| 1309 |
1328 |
_blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
|
| 1310 |
1329 |
left_path.push_back(nca);
|
| 1311 |
1330 |
nca =
|
| 1312 |
1331 |
_blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
|
| 1313 |
1332 |
left_path.push_back(nca);
|
| 1314 |
1333 |
}
|
| 1315 |
1334 |
}
|
| 1316 |
1335 |
}
|
| 1317 |
1336 |
}
|
| 1318 |
1337 |
|
| 1319 |
1338 |
std::vector<int> subblossoms;
|
| 1320 |
1339 |
Arc prev;
|
| 1321 |
1340 |
|
| 1322 |
1341 |
prev = _graph.direct(edge, true);
|
| 1323 |
1342 |
for (int i = 0; left_path[i] != nca; i += 2) {
|
| 1324 |
1343 |
subblossoms.push_back(left_path[i]);
|
| 1325 |
1344 |
(*_blossom_data)[left_path[i]].next = prev;
|
| 1326 |
1345 |
_tree_set->erase(left_path[i]);
|
| 1327 |
1346 |
|
| 1328 |
1347 |
subblossoms.push_back(left_path[i + 1]);
|
| 1329 |
1348 |
(*_blossom_data)[left_path[i + 1]].status = EVEN;
|
| 1330 |
1349 |
oddToEven(left_path[i + 1], tree);
|
| 1331 |
1350 |
_tree_set->erase(left_path[i + 1]);
|
| 1332 |
1351 |
prev = _graph.oppositeArc((*_blossom_data)[left_path[i + 1]].pred);
|
| 1333 |
1352 |
}
|
| 1334 |
1353 |
|
| 1335 |
1354 |
int k = 0;
|
| 1336 |
1355 |
while (right_path[k] != nca) ++k;
|
| 1337 |
1356 |
|
| 1338 |
1357 |
subblossoms.push_back(nca);
|
| 1339 |
1358 |
(*_blossom_data)[nca].next = prev;
|
| 1340 |
1359 |
|
| 1341 |
1360 |
for (int i = k - 2; i >= 0; i -= 2) {
|
| 1342 |
1361 |
subblossoms.push_back(right_path[i + 1]);
|
| 1343 |
1362 |
(*_blossom_data)[right_path[i + 1]].status = EVEN;
|
| 1344 |
1363 |
oddToEven(right_path[i + 1], tree);
|
| 1345 |
1364 |
_tree_set->erase(right_path[i + 1]);
|
| 1346 |
1365 |
|
| 1347 |
1366 |
(*_blossom_data)[right_path[i + 1]].next =
|
| 1348 |
1367 |
(*_blossom_data)[right_path[i + 1]].pred;
|
| 1349 |
1368 |
|
| 1350 |
1369 |
subblossoms.push_back(right_path[i]);
|
| 1351 |
1370 |
_tree_set->erase(right_path[i]);
|
| 1352 |
1371 |
}
|
| 1353 |
1372 |
|
| 1354 |
1373 |
int surface =
|
| 1355 |
1374 |
_blossom_set->join(subblossoms.begin(), subblossoms.end());
|
| 1356 |
1375 |
|
| 1357 |
1376 |
for (int i = 0; i < int(subblossoms.size()); ++i) {
|
| 1358 |
1377 |
if (!_blossom_set->trivial(subblossoms[i])) {
|
| 1359 |
1378 |
(*_blossom_data)[subblossoms[i]].pot += 2 * _delta_sum;
|
| 1360 |
1379 |
}
|
| 1361 |
1380 |
(*_blossom_data)[subblossoms[i]].status = MATCHED;
|
| 1362 |
1381 |
}
|
| 1363 |
1382 |
|
| 1364 |
1383 |
(*_blossom_data)[surface].pot = -2 * _delta_sum;
|
| 1365 |
1384 |
(*_blossom_data)[surface].offset = 0;
|
| 1366 |
1385 |
(*_blossom_data)[surface].status = EVEN;
|
| 1367 |
1386 |
(*_blossom_data)[surface].pred = (*_blossom_data)[nca].pred;
|
| 1368 |
1387 |
(*_blossom_data)[surface].next = (*_blossom_data)[nca].pred;
|
| 1369 |
1388 |
|
| 1370 |
1389 |
_tree_set->insert(surface, tree);
|
| 1371 |
1390 |
_tree_set->erase(nca);
|
| 1372 |
1391 |
}
|
| 1373 |
1392 |
|
| 1374 |
1393 |
void splitBlossom(int blossom) {
|
| 1375 |
1394 |
Arc next = (*_blossom_data)[blossom].next;
|
| 1376 |
1395 |
Arc pred = (*_blossom_data)[blossom].pred;
|
| 1377 |
1396 |
|
| 1378 |
1397 |
int tree = _tree_set->find(blossom);
|
| 1379 |
1398 |
|
| 1380 |
1399 |
(*_blossom_data)[blossom].status = MATCHED;
|
| 1381 |
1400 |
oddToMatched(blossom);
|
| 1382 |
1401 |
if (_delta2->state(blossom) == _delta2->IN_HEAP) {
|
| 1383 |
1402 |
_delta2->erase(blossom);
|
| 1384 |
1403 |
}
|
| 1385 |
1404 |
|
| 1386 |
1405 |
std::vector<int> subblossoms;
|
| 1387 |
1406 |
_blossom_set->split(blossom, std::back_inserter(subblossoms));
|
| 1388 |
1407 |
|
| 1389 |
1408 |
Value offset = (*_blossom_data)[blossom].offset;
|
| 1390 |
1409 |
int b = _blossom_set->find(_graph.source(pred));
|
| 1391 |
1410 |
int d = _blossom_set->find(_graph.source(next));
|
| 1392 |
1411 |
|
| 1393 |
1412 |
int ib = -1, id = -1;
|
| 1394 |
1413 |
for (int i = 0; i < int(subblossoms.size()); ++i) {
|
| 1395 |
1414 |
if (subblossoms[i] == b) ib = i;
|
| 1396 |
1415 |
if (subblossoms[i] == d) id = i;
|
| 1397 |
1416 |
|
| 1398 |
1417 |
(*_blossom_data)[subblossoms[i]].offset = offset;
|
| 1399 |
1418 |
if (!_blossom_set->trivial(subblossoms[i])) {
|
| 1400 |
1419 |
(*_blossom_data)[subblossoms[i]].pot -= 2 * offset;
|
| 1401 |
1420 |
}
|
| 1402 |
1421 |
if (_blossom_set->classPrio(subblossoms[i]) !=
|
| 1403 |
1422 |
std::numeric_limits<Value>::max()) {
|
| 1404 |
1423 |
_delta2->push(subblossoms[i],
|
| 1405 |
1424 |
_blossom_set->classPrio(subblossoms[i]) -
|
| 1406 |
1425 |
(*_blossom_data)[subblossoms[i]].offset);
|
| 1407 |
1426 |
}
|
| 1408 |
1427 |
}
|
| 1409 |
1428 |
|
| 1410 |
1429 |
if (id > ib ? ((id - ib) % 2 == 0) : ((ib - id) % 2 == 1)) {
|
| 1411 |
1430 |
for (int i = (id + 1) % subblossoms.size();
|
| 1412 |
1431 |
i != ib; i = (i + 2) % subblossoms.size()) {
|
| 1413 |
1432 |
int sb = subblossoms[i];
|
| 1414 |
1433 |
int tb = subblossoms[(i + 1) % subblossoms.size()];
|
| 1415 |
1434 |
(*_blossom_data)[sb].next =
|
| 1416 |
1435 |
_graph.oppositeArc((*_blossom_data)[tb].next);
|
| 1417 |
1436 |
}
|
| 1418 |
1437 |
|
| 1419 |
1438 |
for (int i = ib; i != id; i = (i + 2) % subblossoms.size()) {
|
| 1420 |
1439 |
int sb = subblossoms[i];
|
| 1421 |
1440 |
int tb = subblossoms[(i + 1) % subblossoms.size()];
|
| 1422 |
1441 |
int ub = subblossoms[(i + 2) % subblossoms.size()];
|
| 1423 |
1442 |
|
| 1424 |
1443 |
(*_blossom_data)[sb].status = ODD;
|
| 1425 |
1444 |
matchedToOdd(sb);
|
| 1426 |
1445 |
_tree_set->insert(sb, tree);
|
| 1427 |
1446 |
(*_blossom_data)[sb].pred = pred;
|
| 1428 |
1447 |
(*_blossom_data)[sb].next =
|
| 1429 |
1448 |
_graph.oppositeArc((*_blossom_data)[tb].next);
|
| 1430 |
1449 |
|
| 1431 |
1450 |
pred = (*_blossom_data)[ub].next;
|
| 1432 |
1451 |
|
| 1433 |
1452 |
(*_blossom_data)[tb].status = EVEN;
|
| 1434 |
1453 |
matchedToEven(tb, tree);
|
| 1435 |
1454 |
_tree_set->insert(tb, tree);
|
| 1436 |
1455 |
(*_blossom_data)[tb].pred = (*_blossom_data)[tb].next;
|
| 1437 |
1456 |
}
|
| 1438 |
1457 |
|
| 1439 |
1458 |
(*_blossom_data)[subblossoms[id]].status = ODD;
|
| 1440 |
1459 |
matchedToOdd(subblossoms[id]);
|
| 1441 |
1460 |
_tree_set->insert(subblossoms[id], tree);
|
| 1442 |
1461 |
(*_blossom_data)[subblossoms[id]].next = next;
|
| 1443 |
1462 |
(*_blossom_data)[subblossoms[id]].pred = pred;
|
| 1444 |
1463 |
|
| 1445 |
1464 |
} else {
|
| 1446 |
1465 |
|
| 1447 |
1466 |
for (int i = (ib + 1) % subblossoms.size();
|
| 1448 |
1467 |
i != id; i = (i + 2) % subblossoms.size()) {
|
| 1449 |
1468 |
int sb = subblossoms[i];
|
| 1450 |
1469 |
int tb = subblossoms[(i + 1) % subblossoms.size()];
|
| 1451 |
1470 |
(*_blossom_data)[sb].next =
|
| 1452 |
1471 |
_graph.oppositeArc((*_blossom_data)[tb].next);
|
| 1453 |
1472 |
}
|
| 1454 |
1473 |
|
| 1455 |
1474 |
for (int i = id; i != ib; i = (i + 2) % subblossoms.size()) {
|
| 1456 |
1475 |
int sb = subblossoms[i];
|
| 1457 |
1476 |
int tb = subblossoms[(i + 1) % subblossoms.size()];
|
| 1458 |
1477 |
int ub = subblossoms[(i + 2) % subblossoms.size()];
|
| 1459 |
1478 |
|
| 1460 |
1479 |
(*_blossom_data)[sb].status = ODD;
|
| 1461 |
1480 |
matchedToOdd(sb);
|
| 1462 |
1481 |
_tree_set->insert(sb, tree);
|
| 1463 |
1482 |
(*_blossom_data)[sb].next = next;
|
| 1464 |
1483 |
(*_blossom_data)[sb].pred =
|
| 1465 |
1484 |
_graph.oppositeArc((*_blossom_data)[tb].next);
|
| 1466 |
1485 |
|
| 1467 |
1486 |
(*_blossom_data)[tb].status = EVEN;
|
| 1468 |
1487 |
matchedToEven(tb, tree);
|
| 1469 |
1488 |
_tree_set->insert(tb, tree);
|
| 1470 |
1489 |
(*_blossom_data)[tb].pred =
|
| 1471 |
1490 |
(*_blossom_data)[tb].next =
|
| 1472 |
1491 |
_graph.oppositeArc((*_blossom_data)[ub].next);
|
| 1473 |
1492 |
next = (*_blossom_data)[ub].next;
|
| 1474 |
1493 |
}
|
| 1475 |
1494 |
|
| 1476 |
1495 |
(*_blossom_data)[subblossoms[ib]].status = ODD;
|
| 1477 |
1496 |
matchedToOdd(subblossoms[ib]);
|
| 1478 |
1497 |
_tree_set->insert(subblossoms[ib], tree);
|
| 1479 |
1498 |
(*_blossom_data)[subblossoms[ib]].next = next;
|
| 1480 |
1499 |
(*_blossom_data)[subblossoms[ib]].pred = pred;
|
| 1481 |
1500 |
}
|
| 1482 |
1501 |
_tree_set->erase(blossom);
|
| 1483 |
1502 |
}
|
| 1484 |
1503 |
|
| 1485 |
1504 |
void extractBlossom(int blossom, const Node& base, const Arc& matching) {
|
| 1486 |
1505 |
if (_blossom_set->trivial(blossom)) {
|
| 1487 |
1506 |
int bi = (*_node_index)[base];
|
| 1488 |
1507 |
Value pot = (*_node_data)[bi].pot;
|
| 1489 |
1508 |
|
| 1490 |
1509 |
(*_matching)[base] = matching;
|
| 1491 |
1510 |
_blossom_node_list.push_back(base);
|
| 1492 |
1511 |
(*_node_potential)[base] = pot;
|
| 1493 |
1512 |
} else {
|
| 1494 |
1513 |
|
| 1495 |
1514 |
Value pot = (*_blossom_data)[blossom].pot;
|
| 1496 |
1515 |
int bn = _blossom_node_list.size();
|
| 1497 |
1516 |
|
| 1498 |
1517 |
std::vector<int> subblossoms;
|
| 1499 |
1518 |
_blossom_set->split(blossom, std::back_inserter(subblossoms));
|
| 1500 |
1519 |
int b = _blossom_set->find(base);
|
| 1501 |
1520 |
int ib = -1;
|
| 1502 |
1521 |
for (int i = 0; i < int(subblossoms.size()); ++i) {
|
| 1503 |
1522 |
if (subblossoms[i] == b) { ib = i; break; }
|
| 1504 |
1523 |
}
|
| 1505 |
1524 |
|
| 1506 |
1525 |
for (int i = 1; i < int(subblossoms.size()); i += 2) {
|
| 1507 |
1526 |
int sb = subblossoms[(ib + i) % subblossoms.size()];
|
| 1508 |
1527 |
int tb = subblossoms[(ib + i + 1) % subblossoms.size()];
|
| 1509 |
1528 |
|
| 1510 |
1529 |
Arc m = (*_blossom_data)[tb].next;
|
| 1511 |
1530 |
extractBlossom(sb, _graph.target(m), _graph.oppositeArc(m));
|
| 1512 |
1531 |
extractBlossom(tb, _graph.source(m), m);
|
| 1513 |
1532 |
}
|
| 1514 |
1533 |
extractBlossom(subblossoms[ib], base, matching);
|
| 1515 |
1534 |
|
| 1516 |
1535 |
int en = _blossom_node_list.size();
|
| 1517 |
1536 |
|
| 1518 |
1537 |
_blossom_potential.push_back(BlossomVariable(bn, en, pot));
|
| 1519 |
1538 |
}
|
| 1520 |
1539 |
}
|
| 1521 |
1540 |
|
| 1522 |
1541 |
void extractMatching() {
|
| 1523 |
1542 |
std::vector<int> blossoms;
|
| 1524 |
1543 |
for (typename BlossomSet::ClassIt c(*_blossom_set); c != INVALID; ++c) {
|
| 1525 |
1544 |
blossoms.push_back(c);
|
| 1526 |
1545 |
}
|
| 1527 |
1546 |
|
| 1528 |
1547 |
for (int i = 0; i < int(blossoms.size()); ++i) {
|
| 1529 |
1548 |
if ((*_blossom_data)[blossoms[i]].next != INVALID) {
|
| 1530 |
1549 |
|
| 1531 |
1550 |
Value offset = (*_blossom_data)[blossoms[i]].offset;
|
| 1532 |
1551 |
(*_blossom_data)[blossoms[i]].pot += 2 * offset;
|
| 1533 |
1552 |
for (typename BlossomSet::ItemIt n(*_blossom_set, blossoms[i]);
|
| 1534 |
1553 |
n != INVALID; ++n) {
|
| 1535 |
1554 |
(*_node_data)[(*_node_index)[n]].pot -= offset;
|
| 1536 |
1555 |
}
|
| 1537 |
1556 |
|
| 1538 |
1557 |
Arc matching = (*_blossom_data)[blossoms[i]].next;
|
| 1539 |
1558 |
Node base = _graph.source(matching);
|
| 1540 |
1559 |
extractBlossom(blossoms[i], base, matching);
|
| 1541 |
1560 |
} else {
|
| 1542 |
1561 |
Node base = (*_blossom_data)[blossoms[i]].base;
|
| 1543 |
1562 |
extractBlossom(blossoms[i], base, INVALID);
|
| 1544 |
1563 |
}
|
| 1545 |
1564 |
}
|
| 1546 |
1565 |
}
|
| 1547 |
1566 |
|
| 1548 |
1567 |
public:
|
| 1549 |
1568 |
|
| 1550 |
1569 |
/// \brief Constructor
|
| 1551 |
1570 |
///
|
| 1552 |
1571 |
/// Constructor.
|
| 1553 |
1572 |
MaxWeightedMatching(const Graph& graph, const WeightMap& weight)
|
| 1554 |
1573 |
: _graph(graph), _weight(weight), _matching(0),
|
| 1555 |
1574 |
_node_potential(0), _blossom_potential(), _blossom_node_list(),
|
| 1556 |
1575 |
_node_num(0), _blossom_num(0),
|
| 1557 |
1576 |
|
| 1558 |
1577 |
_blossom_index(0), _blossom_set(0), _blossom_data(0),
|
| 1559 |
1578 |
_node_index(0), _node_heap_index(0), _node_data(0),
|
| 1560 |
1579 |
_tree_set_index(0), _tree_set(0),
|
| 1561 |
1580 |
|
| 1562 |
1581 |
_delta1_index(0), _delta1(0),
|
| 1563 |
1582 |
_delta2_index(0), _delta2(0),
|
| 1564 |
1583 |
_delta3_index(0), _delta3(0),
|
| 1565 |
1584 |
_delta4_index(0), _delta4(0),
|
| 1566 |
1585 |
|
| 1567 |
1586 |
_delta_sum(), _unmatched(0),
|
| 1568 |
1587 |
|
| 1569 |
1588 |
_fractional(0)
|
| 1570 |
1589 |
{}
|
| 1571 |
1590 |
|
| 1572 |
1591 |
~MaxWeightedMatching() {
|
| 1573 |
1592 |
destroyStructures();
|
| 1574 |
1593 |
if (_fractional) {
|
| 1575 |
1594 |
delete _fractional;
|
| 1576 |
1595 |
}
|
| 1577 |
1596 |
}
|
| 1578 |
1597 |
|
| 1579 |
1598 |
/// \name Execution Control
|
| 1580 |
1599 |
/// The simplest way to execute the algorithm is to use the
|
| 1581 |
1600 |
/// \ref run() member function.
|
| 1582 |
1601 |
|
| 1583 |
1602 |
///@{
|
| 1584 |
1603 |
|
| 1585 |
1604 |
/// \brief Initialize the algorithm
|
| 1586 |
1605 |
///
|
| 1587 |
1606 |
/// This function initializes the algorithm.
|
| 1588 |
1607 |
void init() {
|
| 1589 |
1608 |
createStructures();
|
| 1590 |
1609 |
|
|
1610 |
_blossom_node_list.clear();
|
|
1611 |
_blossom_potential.clear();
|
|
1612 |
|
| 1591 |
1613 |
for (ArcIt e(_graph); e != INVALID; ++e) {
|
| 1592 |
1614 |
(*_node_heap_index)[e] = BinHeap<Value, IntArcMap>::PRE_HEAP;
|
| 1593 |
1615 |
}
|
| 1594 |
1616 |
for (NodeIt n(_graph); n != INVALID; ++n) {
|
| 1595 |
1617 |
(*_delta1_index)[n] = _delta1->PRE_HEAP;
|
| 1596 |
1618 |
}
|
| 1597 |
1619 |
for (EdgeIt e(_graph); e != INVALID; ++e) {
|
| 1598 |
1620 |
(*_delta3_index)[e] = _delta3->PRE_HEAP;
|
| 1599 |
1621 |
}
|
| 1600 |
1622 |
for (int i = 0; i < _blossom_num; ++i) {
|
| 1601 |
1623 |
(*_delta2_index)[i] = _delta2->PRE_HEAP;
|
| 1602 |
1624 |
(*_delta4_index)[i] = _delta4->PRE_HEAP;
|
| 1603 |
1625 |
}
|
| 1604 |
|
|
|
1626 |
|
| 1605 |
1627 |
_unmatched = _node_num;
|
| 1606 |
1628 |
|
|
1629 |
_delta1->clear();
|
|
1630 |
_delta2->clear();
|
|
1631 |
_delta3->clear();
|
|
1632 |
_delta4->clear();
|
|
1633 |
_blossom_set->clear();
|
|
1634 |
_tree_set->clear();
|
|
1635 |
|
| 1607 |
1636 |
int index = 0;
|
| 1608 |
1637 |
for (NodeIt n(_graph); n != INVALID; ++n) {
|
| 1609 |
1638 |
Value max = 0;
|
| 1610 |
1639 |
for (OutArcIt e(_graph, n); e != INVALID; ++e) {
|
| 1611 |
1640 |
if (_graph.target(e) == n) continue;
|
| 1612 |
1641 |
if ((dualScale * _weight[e]) / 2 > max) {
|
| 1613 |
1642 |
max = (dualScale * _weight[e]) / 2;
|
| 1614 |
1643 |
}
|
| 1615 |
1644 |
}
|
| 1616 |
1645 |
(*_node_index)[n] = index;
|
|
1646 |
(*_node_data)[index].heap_index.clear();
|
|
1647 |
(*_node_data)[index].heap.clear();
|
| 1617 |
1648 |
(*_node_data)[index].pot = max;
|
| 1618 |
1649 |
_delta1->push(n, max);
|
| 1619 |
1650 |
int blossom =
|
| 1620 |
1651 |
_blossom_set->insert(n, std::numeric_limits<Value>::max());
|
| 1621 |
1652 |
|
| 1622 |
1653 |
_tree_set->insert(blossom);
|
| 1623 |
1654 |
|
| 1624 |
1655 |
(*_blossom_data)[blossom].status = EVEN;
|
| 1625 |
1656 |
(*_blossom_data)[blossom].pred = INVALID;
|
| 1626 |
1657 |
(*_blossom_data)[blossom].next = INVALID;
|
| 1627 |
1658 |
(*_blossom_data)[blossom].pot = 0;
|
| 1628 |
1659 |
(*_blossom_data)[blossom].offset = 0;
|
| 1629 |
1660 |
++index;
|
| 1630 |
1661 |
}
|
| 1631 |
1662 |
for (EdgeIt e(_graph); e != INVALID; ++e) {
|
| 1632 |
1663 |
int si = (*_node_index)[_graph.u(e)];
|
| 1633 |
1664 |
int ti = (*_node_index)[_graph.v(e)];
|
| 1634 |
1665 |
if (_graph.u(e) != _graph.v(e)) {
|
| 1635 |
1666 |
_delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot -
|
| 1636 |
1667 |
dualScale * _weight[e]) / 2);
|
| 1637 |
1668 |
}
|
| 1638 |
1669 |
}
|
| 1639 |
1670 |
}
|
| 1640 |
1671 |
|
| 1641 |
1672 |
/// \brief Initialize the algorithm with fractional matching
|
| 1642 |
1673 |
///
|
| 1643 |
1674 |
/// This function initializes the algorithm with a fractional
|
| 1644 |
1675 |
/// matching. This initialization is also called jumpstart heuristic.
|
| 1645 |
1676 |
void fractionalInit() {
|
| 1646 |
1677 |
createStructures();
|
| 1647 |
1678 |
|
| 1648 |
1679 |
if (_fractional == 0) {
|
| 1649 |
1680 |
_fractional = new FractionalMatching(_graph, _weight, false);
|
| 1650 |
1681 |
}
|
| 1651 |
1682 |
_fractional->run();
|
| 1652 |
1683 |
|
| 1653 |
1684 |
for (ArcIt e(_graph); e != INVALID; ++e) {
|
| 1654 |
1685 |
(*_node_heap_index)[e] = BinHeap<Value, IntArcMap>::PRE_HEAP;
|
| 1655 |
1686 |
}
|
| 1656 |
1687 |
for (NodeIt n(_graph); n != INVALID; ++n) {
|
| 1657 |
1688 |
(*_delta1_index)[n] = _delta1->PRE_HEAP;
|
| 1658 |
1689 |
}
|
| 1659 |
1690 |
for (EdgeIt e(_graph); e != INVALID; ++e) {
|
| 1660 |
1691 |
(*_delta3_index)[e] = _delta3->PRE_HEAP;
|
| 1661 |
1692 |
}
|
| 1662 |
1693 |
for (int i = 0; i < _blossom_num; ++i) {
|
| 1663 |
1694 |
(*_delta2_index)[i] = _delta2->PRE_HEAP;
|
| 1664 |
1695 |
(*_delta4_index)[i] = _delta4->PRE_HEAP;
|
| 1665 |
1696 |
}
|
| 1666 |
1697 |
|
| 1667 |
1698 |
_unmatched = 0;
|
| 1668 |
1699 |
|
| 1669 |
1700 |
int index = 0;
|
| 1670 |
1701 |
for (NodeIt n(_graph); n != INVALID; ++n) {
|
| 1671 |
1702 |
Value pot = _fractional->nodeValue(n);
|
| 1672 |
1703 |
(*_node_index)[n] = index;
|
| 1673 |
1704 |
(*_node_data)[index].pot = pot;
|
| 1674 |
1705 |
int blossom =
|
| 1675 |
1706 |
_blossom_set->insert(n, std::numeric_limits<Value>::max());
|
| 1676 |
1707 |
|
| 1677 |
1708 |
(*_blossom_data)[blossom].status = MATCHED;
|
| 1678 |
1709 |
(*_blossom_data)[blossom].pred = INVALID;
|
| 1679 |
1710 |
(*_blossom_data)[blossom].next = _fractional->matching(n);
|
| 1680 |
1711 |
if (_fractional->matching(n) == INVALID) {
|
| 1681 |
1712 |
(*_blossom_data)[blossom].base = n;
|
| 1682 |
1713 |
}
|
| 1683 |
1714 |
(*_blossom_data)[blossom].pot = 0;
|
| 1684 |
1715 |
(*_blossom_data)[blossom].offset = 0;
|
| 1685 |
1716 |
++index;
|
| 1686 |
1717 |
}
|
| 1687 |
1718 |
|
| 1688 |
1719 |
typename Graph::template NodeMap<bool> processed(_graph, false);
|
| 1689 |
1720 |
for (NodeIt n(_graph); n != INVALID; ++n) {
|
| 1690 |
1721 |
if (processed[n]) continue;
|
| 1691 |
1722 |
processed[n] = true;
|
| 1692 |
1723 |
if (_fractional->matching(n) == INVALID) continue;
|
| 1693 |
1724 |
int num = 1;
|
| 1694 |
1725 |
Node v = _graph.target(_fractional->matching(n));
|
| 1695 |
1726 |
while (n != v) {
|
| 1696 |
1727 |
processed[v] = true;
|
| 1697 |
1728 |
v = _graph.target(_fractional->matching(v));
|
| 1698 |
1729 |
++num;
|
| 1699 |
1730 |
}
|
| 1700 |
1731 |
|
| 1701 |
1732 |
if (num % 2 == 1) {
|
| 1702 |
1733 |
std::vector<int> subblossoms(num);
|
| 1703 |
1734 |
|
| 1704 |
1735 |
subblossoms[--num] = _blossom_set->find(n);
|
| 1705 |
1736 |
_delta1->push(n, _fractional->nodeValue(n));
|
| 1706 |
1737 |
v = _graph.target(_fractional->matching(n));
|
| 1707 |
1738 |
while (n != v) {
|
| 1708 |
1739 |
subblossoms[--num] = _blossom_set->find(v);
|
| 1709 |
1740 |
_delta1->push(v, _fractional->nodeValue(v));
|
| 1710 |
1741 |
v = _graph.target(_fractional->matching(v));
|
| 1711 |
1742 |
}
|
| 1712 |
1743 |
|
| 1713 |
1744 |
int surface =
|
| 1714 |
1745 |
_blossom_set->join(subblossoms.begin(), subblossoms.end());
|
| 1715 |
1746 |
(*_blossom_data)[surface].status = EVEN;
|
| 1716 |
1747 |
(*_blossom_data)[surface].pred = INVALID;
|
| 1717 |
1748 |
(*_blossom_data)[surface].next = INVALID;
|
| 1718 |
1749 |
(*_blossom_data)[surface].pot = 0;
|
| 1719 |
1750 |
(*_blossom_data)[surface].offset = 0;
|
| 1720 |
1751 |
|
| 1721 |
1752 |
_tree_set->insert(surface);
|
| 1722 |
1753 |
++_unmatched;
|
| 1723 |
1754 |
}
|
| 1724 |
1755 |
}
|
| 1725 |
1756 |
|
| 1726 |
1757 |
for (EdgeIt e(_graph); e != INVALID; ++e) {
|
| 1727 |
1758 |
int si = (*_node_index)[_graph.u(e)];
|
| 1728 |
1759 |
int sb = _blossom_set->find(_graph.u(e));
|
| 1729 |
1760 |
int ti = (*_node_index)[_graph.v(e)];
|
| 1730 |
1761 |
int tb = _blossom_set->find(_graph.v(e));
|
| 1731 |
1762 |
if ((*_blossom_data)[sb].status == EVEN &&
|
| 1732 |
1763 |
(*_blossom_data)[tb].status == EVEN && sb != tb) {
|
| 1733 |
1764 |
_delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot -
|
| 1734 |
1765 |
dualScale * _weight[e]) / 2);
|
| 1735 |
1766 |
}
|
| 1736 |
1767 |
}
|
| 1737 |
1768 |
|
| 1738 |
1769 |
for (NodeIt n(_graph); n != INVALID; ++n) {
|
| 1739 |
1770 |
int nb = _blossom_set->find(n);
|
| 1740 |
1771 |
if ((*_blossom_data)[nb].status != MATCHED) continue;
|
| 1741 |
1772 |
int ni = (*_node_index)[n];
|
| 1742 |
1773 |
|
| 1743 |
1774 |
for (OutArcIt e(_graph, n); e != INVALID; ++e) {
|
| 1744 |
1775 |
Node v = _graph.target(e);
|
| 1745 |
1776 |
int vb = _blossom_set->find(v);
|
| 1746 |
1777 |
int vi = (*_node_index)[v];
|
| 1747 |
1778 |
|
| 1748 |
1779 |
Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
|
| 1749 |
1780 |
dualScale * _weight[e];
|
| 1750 |
1781 |
|
| 1751 |
1782 |
if ((*_blossom_data)[vb].status == EVEN) {
|
| 1752 |
1783 |
|
| 1753 |
1784 |
int vt = _tree_set->find(vb);
|
| 1754 |
1785 |
|
| 1755 |
1786 |
typename std::map<int, Arc>::iterator it =
|
| 1756 |
1787 |
(*_node_data)[ni].heap_index.find(vt);
|
| 1757 |
1788 |
|
| 1758 |
1789 |
if (it != (*_node_data)[ni].heap_index.end()) {
|
| 1759 |
1790 |
if ((*_node_data)[ni].heap[it->second] > rw) {
|
| 1760 |
1791 |
(*_node_data)[ni].heap.replace(it->second, e);
|
| 1761 |
1792 |
(*_node_data)[ni].heap.decrease(e, rw);
|
| 1762 |
1793 |
it->second = e;
|
| 1763 |
1794 |
}
|
| 1764 |
1795 |
} else {
|
| 1765 |
1796 |
(*_node_data)[ni].heap.push(e, rw);
|
| 1766 |
1797 |
(*_node_data)[ni].heap_index.insert(std::make_pair(vt, e));
|
| 1767 |
1798 |
}
|
| 1768 |
1799 |
}
|
| 1769 |
1800 |
}
|
| 1770 |
1801 |
|
| 1771 |
1802 |
if (!(*_node_data)[ni].heap.empty()) {
|
| 1772 |
1803 |
_blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
|
| 1773 |
1804 |
_delta2->push(nb, _blossom_set->classPrio(nb));
|
| 1774 |
1805 |
}
|
| 1775 |
1806 |
}
|
| 1776 |
1807 |
}
|
| 1777 |
1808 |
|
| 1778 |
1809 |
/// \brief Start the algorithm
|
| 1779 |
1810 |
///
|
| 1780 |
1811 |
/// This function starts the algorithm.
|
| 1781 |
1812 |
///
|
| 1782 |
1813 |
/// \pre \ref init() or \ref fractionalInit() must be called
|
| 1783 |
1814 |
/// before using this function.
|
| 1784 |
1815 |
void start() {
|
| 1785 |
1816 |
enum OpType {
|
| 1786 |
1817 |
D1, D2, D3, D4
|
| 1787 |
1818 |
};
|
| 1788 |
1819 |
|
| 1789 |
1820 |
while (_unmatched > 0) {
|
| 1790 |
1821 |
Value d1 = !_delta1->empty() ?
|
| 1791 |
1822 |
_delta1->prio() : std::numeric_limits<Value>::max();
|
| 1792 |
1823 |
|
| 1793 |
1824 |
Value d2 = !_delta2->empty() ?
|
| 1794 |
1825 |
_delta2->prio() : std::numeric_limits<Value>::max();
|
| 1795 |
1826 |
|
| 1796 |
1827 |
Value d3 = !_delta3->empty() ?
|
| 1797 |
1828 |
_delta3->prio() : std::numeric_limits<Value>::max();
|
| 1798 |
1829 |
|
| 1799 |
1830 |
Value d4 = !_delta4->empty() ?
|
| 1800 |
1831 |
_delta4->prio() : std::numeric_limits<Value>::max();
|
| 1801 |
1832 |
|
| 1802 |
1833 |
_delta_sum = d3; OpType ot = D3;
|
| 1803 |
1834 |
if (d1 < _delta_sum) { _delta_sum = d1; ot = D1; }
|
| 1804 |
1835 |
if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; }
|
| 1805 |
1836 |
if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; }
|
| 1806 |
1837 |
|
| 1807 |
1838 |
switch (ot) {
|
| 1808 |
1839 |
case D1:
|
| 1809 |
1840 |
{
|
| 1810 |
1841 |
Node n = _delta1->top();
|
| 1811 |
1842 |
unmatchNode(n);
|
| 1812 |
1843 |
--_unmatched;
|
| 1813 |
1844 |
}
|
| 1814 |
1845 |
break;
|
| 1815 |
1846 |
case D2:
|
| 1816 |
1847 |
{
|
| 1817 |
1848 |
int blossom = _delta2->top();
|
| 1818 |
1849 |
Node n = _blossom_set->classTop(blossom);
|
| 1819 |
1850 |
Arc a = (*_node_data)[(*_node_index)[n]].heap.top();
|
| 1820 |
1851 |
if ((*_blossom_data)[blossom].next == INVALID) {
|
| 1821 |
1852 |
augmentOnArc(a);
|
| 1822 |
1853 |
--_unmatched;
|
| 1823 |
1854 |
} else {
|
| 1824 |
1855 |
extendOnArc(a);
|
| 1825 |
1856 |
}
|
| 1826 |
1857 |
}
|
| 1827 |
1858 |
break;
|
| 1828 |
1859 |
case D3:
|
| 1829 |
1860 |
{
|
| 1830 |
1861 |
Edge e = _delta3->top();
|
| 1831 |
1862 |
|
| 1832 |
1863 |
int left_blossom = _blossom_set->find(_graph.u(e));
|
| 1833 |
1864 |
int right_blossom = _blossom_set->find(_graph.v(e));
|
| 1834 |
1865 |
|
| 1835 |
1866 |
if (left_blossom == right_blossom) {
|
| 1836 |
1867 |
_delta3->pop();
|
| 1837 |
1868 |
} else {
|
| 1838 |
1869 |
int left_tree = _tree_set->find(left_blossom);
|
| 1839 |
1870 |
int right_tree = _tree_set->find(right_blossom);
|
| 1840 |
1871 |
|
| 1841 |
1872 |
if (left_tree == right_tree) {
|
| 1842 |
1873 |
shrinkOnEdge(e, left_tree);
|
| 1843 |
1874 |
} else {
|
| 1844 |
1875 |
augmentOnEdge(e);
|
| 1845 |
1876 |
_unmatched -= 2;
|
| 1846 |
1877 |
}
|
| 1847 |
1878 |
}
|
| 1848 |
1879 |
} break;
|
| 1849 |
1880 |
case D4:
|
| 1850 |
1881 |
splitBlossom(_delta4->top());
|
| 1851 |
1882 |
break;
|
| 1852 |
1883 |
}
|
| 1853 |
1884 |
}
|
| 1854 |
1885 |
extractMatching();
|
| 1855 |
1886 |
}
|
| 1856 |
1887 |
|
| 1857 |
1888 |
/// \brief Run the algorithm.
|
| 1858 |
1889 |
///
|
| 1859 |
1890 |
/// This method runs the \c %MaxWeightedMatching algorithm.
|
| 1860 |
1891 |
///
|
| 1861 |
1892 |
/// \note mwm.run() is just a shortcut of the following code.
|
| 1862 |
1893 |
/// \code
|
| 1863 |
1894 |
/// mwm.fractionalInit();
|
| 1864 |
1895 |
/// mwm.start();
|
| 1865 |
1896 |
/// \endcode
|
| 1866 |
1897 |
void run() {
|
| 1867 |
1898 |
fractionalInit();
|
| 1868 |
1899 |
start();
|
| 1869 |
1900 |
}
|
| 1870 |
1901 |
|
| 1871 |
1902 |
/// @}
|
| 1872 |
1903 |
|
| 1873 |
1904 |
/// \name Primal Solution
|
| 1874 |
1905 |
/// Functions to get the primal solution, i.e. the maximum weighted
|
| 1875 |
1906 |
/// matching.\n
|
| 1876 |
1907 |
/// Either \ref run() or \ref start() function should be called before
|
| 1877 |
1908 |
/// using them.
|
| 1878 |
1909 |
|
| 1879 |
1910 |
/// @{
|
| 1880 |
1911 |
|
| 1881 |
1912 |
/// \brief Return the weight of the matching.
|
| 1882 |
1913 |
///
|
| 1883 |
1914 |
/// This function returns the weight of the found matching.
|
| 1884 |
1915 |
///
|
| 1885 |
1916 |
/// \pre Either run() or start() must be called before using this function.
|
| 1886 |
1917 |
Value matchingWeight() const {
|
| 1887 |
1918 |
Value sum = 0;
|
| 1888 |
1919 |
for (NodeIt n(_graph); n != INVALID; ++n) {
|
| 1889 |
1920 |
if ((*_matching)[n] != INVALID) {
|
| 1890 |
1921 |
sum += _weight[(*_matching)[n]];
|
| 1891 |
1922 |
}
|
| 1892 |
1923 |
}
|
| 1893 |
1924 |
return sum / 2;
|
| 1894 |
1925 |
}
|
| 1895 |
1926 |
|
| 1896 |
1927 |
/// \brief Return the size (cardinality) of the matching.
|
| 1897 |
1928 |
///
|
| 1898 |
1929 |
/// This function returns the size (cardinality) of the found matching.
|
| 1899 |
1930 |
///
|
| 1900 |
1931 |
/// \pre Either run() or start() must be called before using this function.
|
| 1901 |
1932 |
int matchingSize() const {
|
| 1902 |
1933 |
int num = 0;
|
| 1903 |
1934 |
for (NodeIt n(_graph); n != INVALID; ++n) {
|
| 1904 |
1935 |
if ((*_matching)[n] != INVALID) {
|
| 1905 |
1936 |
++num;
|
| 1906 |
1937 |
}
|
| 1907 |
1938 |
}
|
| 1908 |
1939 |
return num /= 2;
|
| 1909 |
1940 |
}
|
| 1910 |
1941 |
|
| 1911 |
1942 |
/// \brief Return \c true if the given edge is in the matching.
|
| 1912 |
1943 |
///
|
| 1913 |
1944 |
/// This function returns \c true if the given edge is in the found
|
| 1914 |
1945 |
/// matching.
|
| 1915 |
1946 |
///
|
| 1916 |
1947 |
/// \pre Either run() or start() must be called before using this function.
|
| 1917 |
1948 |
bool matching(const Edge& edge) const {
|
| 1918 |
1949 |
return edge == (*_matching)[_graph.u(edge)];
|
| 1919 |
1950 |
}
|
| 1920 |
1951 |
|
| 1921 |
1952 |
/// \brief Return the matching arc (or edge) incident to the given node.
|
| 1922 |
1953 |
///
|
| 1923 |
1954 |
/// This function returns the matching arc (or edge) incident to the
|
| 1924 |
1955 |
/// given node in the found matching or \c INVALID if the node is
|
| 1925 |
1956 |
/// not covered by the matching.
|
| 1926 |
1957 |
///
|
| 1927 |
1958 |
/// \pre Either run() or start() must be called before using this function.
|
| 1928 |
1959 |
Arc matching(const Node& node) const {
|
| 1929 |
1960 |
return (*_matching)[node];
|
| 1930 |
1961 |
}
|
| 1931 |
1962 |
|
| 1932 |
1963 |
/// \brief Return a const reference to the matching map.
|
| 1933 |
1964 |
///
|
| 1934 |
1965 |
/// This function returns a const reference to a node map that stores
|
| 1935 |
1966 |
/// the matching arc (or edge) incident to each node.
|
| 1936 |
1967 |
const MatchingMap& matchingMap() const {
|
| 1937 |
1968 |
return *_matching;
|
| 1938 |
1969 |
}
|
| 1939 |
1970 |
|
| 1940 |
1971 |
/// \brief Return the mate of the given node.
|
| 1941 |
1972 |
///
|
| 1942 |
1973 |
/// This function returns the mate of the given node in the found
|
| 1943 |
1974 |
/// matching or \c INVALID if the node is not covered by the matching.
|
| 1944 |
1975 |
///
|
| 1945 |
1976 |
/// \pre Either run() or start() must be called before using this function.
|
| 1946 |
1977 |
Node mate(const Node& node) const {
|
| 1947 |
1978 |
return (*_matching)[node] != INVALID ?
|
| 1948 |
1979 |
_graph.target((*_matching)[node]) : INVALID;
|
| 1949 |
1980 |
}
|
| 1950 |
1981 |
|
| 1951 |
1982 |
/// @}
|
| 1952 |
1983 |
|
| 1953 |
1984 |
/// \name Dual Solution
|
| 1954 |
1985 |
/// Functions to get the dual solution.\n
|
| 1955 |
1986 |
/// Either \ref run() or \ref start() function should be called before
|
| 1956 |
1987 |
/// using them.
|
| 1957 |
1988 |
|
| 1958 |
1989 |
/// @{
|
| 1959 |
1990 |
|
| 1960 |
1991 |
/// \brief Return the value of the dual solution.
|
| 1961 |
1992 |
///
|
| 1962 |
1993 |
/// This function returns the value of the dual solution.
|
| 1963 |
1994 |
/// It should be equal to the primal value scaled by \ref dualScale
|
| 1964 |
1995 |
/// "dual scale".
|
| 1965 |
1996 |
///
|
| 1966 |
1997 |
/// \pre Either run() or start() must be called before using this function.
|
| 1967 |
1998 |
Value dualValue() const {
|
| 1968 |
1999 |
Value sum = 0;
|
| 1969 |
2000 |
for (NodeIt n(_graph); n != INVALID; ++n) {
|
| 1970 |
2001 |
sum += nodeValue(n);
|
| 1971 |
2002 |
}
|
| 1972 |
2003 |
for (int i = 0; i < blossomNum(); ++i) {
|
| 1973 |
2004 |
sum += blossomValue(i) * (blossomSize(i) / 2);
|
| 1974 |
2005 |
}
|
| 1975 |
2006 |
return sum;
|
| 1976 |
2007 |
}
|
| 1977 |
2008 |
|
| 1978 |
2009 |
/// \brief Return the dual value (potential) of the given node.
|
| 1979 |
2010 |
///
|
| 1980 |
2011 |
/// This function returns the dual value (potential) of the given node.
|
| 1981 |
2012 |
///
|
| 1982 |
2013 |
/// \pre Either run() or start() must be called before using this function.
|
| 1983 |
2014 |
Value nodeValue(const Node& n) const {
|
| 1984 |
2015 |
return (*_node_potential)[n];
|
| 1985 |
2016 |
}
|
| 1986 |
2017 |
|
| 1987 |
2018 |
/// \brief Return the number of the blossoms in the basis.
|
| 1988 |
2019 |
///
|
| 1989 |
2020 |
/// This function returns the number of the blossoms in the basis.
|
| 1990 |
2021 |
///
|
| 1991 |
2022 |
/// \pre Either run() or start() must be called before using this function.
|
| 1992 |
2023 |
/// \see BlossomIt
|
| 1993 |
2024 |
int blossomNum() const {
|
| 1994 |
2025 |
return _blossom_potential.size();
|
| 1995 |
2026 |
}
|
| 1996 |
2027 |
|
| 1997 |
2028 |
/// \brief Return the number of the nodes in the given blossom.
|
| 1998 |
2029 |
///
|
| 1999 |
2030 |
/// This function returns the number of the nodes in the given blossom.
|
| 2000 |
2031 |
///
|
| 2001 |
2032 |
/// \pre Either run() or start() must be called before using this function.
|
| 2002 |
2033 |
/// \see BlossomIt
|
| 2003 |
2034 |
int blossomSize(int k) const {
|
| 2004 |
2035 |
return _blossom_potential[k].end - _blossom_potential[k].begin;
|
| 2005 |
2036 |
}
|
| 2006 |
2037 |
|
| 2007 |
2038 |
/// \brief Return the dual value (ptential) of the given blossom.
|
| 2008 |
2039 |
///
|
| 2009 |
2040 |
/// This function returns the dual value (ptential) of the given blossom.
|
| 2010 |
2041 |
///
|
| 2011 |
2042 |
/// \pre Either run() or start() must be called before using this function.
|
| 2012 |
2043 |
Value blossomValue(int k) const {
|
| 2013 |
2044 |
return _blossom_potential[k].value;
|
| 2014 |
2045 |
}
|
| 2015 |
2046 |
|
| 2016 |
2047 |
/// \brief Iterator for obtaining the nodes of a blossom.
|
| 2017 |
2048 |
///
|
| 2018 |
2049 |
/// This class provides an iterator for obtaining the nodes of the
|
| 2019 |
2050 |
/// given blossom. It lists a subset of the nodes.
|
| 2020 |
2051 |
/// Before using this iterator, you must allocate a
|
| 2021 |
2052 |
/// MaxWeightedMatching class and execute it.
|
| 2022 |
2053 |
class BlossomIt {
|
| 2023 |
2054 |
public:
|
| 2024 |
2055 |
|
| 2025 |
2056 |
/// \brief Constructor.
|
| 2026 |
2057 |
///
|
| 2027 |
2058 |
/// Constructor to get the nodes of the given variable.
|
| 2028 |
2059 |
///
|
| 2029 |
2060 |
/// \pre Either \ref MaxWeightedMatching::run() "algorithm.run()" or
|
| 2030 |
2061 |
/// \ref MaxWeightedMatching::start() "algorithm.start()" must be
|
| 2031 |
2062 |
/// called before initializing this iterator.
|
| 2032 |
2063 |
BlossomIt(const MaxWeightedMatching& algorithm, int variable)
|
| 2033 |
2064 |
: _algorithm(&algorithm)
|
| 2034 |
2065 |
{
|
| 2035 |
2066 |
_index = _algorithm->_blossom_potential[variable].begin;
|
| 2036 |
2067 |
_last = _algorithm->_blossom_potential[variable].end;
|
| 2037 |
2068 |
}
|
| 2038 |
2069 |
|
| 2039 |
2070 |
/// \brief Conversion to \c Node.
|
| 2040 |
2071 |
///
|
| 2041 |
2072 |
/// Conversion to \c Node.
|
| 2042 |
2073 |
operator Node() const {
|
| 2043 |
2074 |
return _algorithm->_blossom_node_list[_index];
|
| 2044 |
2075 |
}
|
| 2045 |
2076 |
|
| 2046 |
2077 |
/// \brief Increment operator.
|
| 2047 |
2078 |
///
|
| 2048 |
2079 |
/// Increment operator.
|
| 2049 |
2080 |
BlossomIt& operator++() {
|
| 2050 |
2081 |
++_index;
|
| 2051 |
2082 |
return *this;
|
| 2052 |
2083 |
}
|
| 2053 |
2084 |
|
| 2054 |
2085 |
/// \brief Validity checking
|
| 2055 |
2086 |
///
|
| 2056 |
2087 |
/// Checks whether the iterator is invalid.
|
| 2057 |
2088 |
bool operator==(Invalid) const { return _index == _last; }
|
| 2058 |
2089 |
|
| 2059 |
2090 |
/// \brief Validity checking
|
| 2060 |
2091 |
///
|
| 2061 |
2092 |
/// Checks whether the iterator is valid.
|
| 2062 |
2093 |
bool operator!=(Invalid) const { return _index != _last; }
|
| 2063 |
2094 |
|
| 2064 |
2095 |
private:
|
| 2065 |
2096 |
const MaxWeightedMatching* _algorithm;
|
| 2066 |
2097 |
int _last;
|
| 2067 |
2098 |
int _index;
|
| 2068 |
2099 |
};
|
| 2069 |
2100 |
|
| 2070 |
2101 |
/// @}
|
| 2071 |
2102 |
|
| 2072 |
2103 |
};
|
| 2073 |
2104 |
|
| 2074 |
2105 |
/// \ingroup matching
|
| 2075 |
2106 |
///
|
| 2076 |
2107 |
/// \brief Weighted perfect matching in general graphs
|
| 2077 |
2108 |
///
|
| 2078 |
2109 |
/// This class provides an efficient implementation of Edmond's
|
| 2079 |
2110 |
/// maximum weighted perfect matching algorithm. The implementation
|
| 2080 |
2111 |
/// is based on extensive use of priority queues and provides
|
| 2081 |
2112 |
/// \f$O(nm\log n)\f$ time complexity.
|
| 2082 |
2113 |
///
|
| 2083 |
2114 |
/// The maximum weighted perfect matching problem is to find a subset of
|
| 2084 |
2115 |
/// the edges in an undirected graph with maximum overall weight for which
|
| 2085 |
2116 |
/// each node has exactly one incident edge.
|
| 2086 |
2117 |
/// It can be formulated with the following linear program.
|
| 2087 |
2118 |
/// \f[ \sum_{e \in \delta(u)}x_e = 1 \quad \forall u\in V\f]
|
| 2088 |
2119 |
/** \f[ \sum_{e \in \gamma(B)}x_e \le \frac{\vert B \vert - 1}{2}
|
| 2089 |
2120 |
\quad \forall B\in\mathcal{O}\f] */
|
| 2090 |
2121 |
/// \f[x_e \ge 0\quad \forall e\in E\f]
|
| 2091 |
2122 |
/// \f[\max \sum_{e\in E}x_ew_e\f]
|
| 2092 |
2123 |
/// where \f$\delta(X)\f$ is the set of edges incident to a node in
|
| 2093 |
2124 |
/// \f$X\f$, \f$\gamma(X)\f$ is the set of edges with both ends in
|
| 2094 |
2125 |
/// \f$X\f$ and \f$\mathcal{O}\f$ is the set of odd cardinality
|
| 2095 |
2126 |
/// subsets of the nodes.
|
| 2096 |
2127 |
///
|
| 2097 |
2128 |
/// The algorithm calculates an optimal matching and a proof of the
|
| 2098 |
2129 |
/// optimality. The solution of the dual problem can be used to check
|
| 2099 |
2130 |
/// the result of the algorithm. The dual linear problem is the
|
| 2100 |
2131 |
/// following.
|
| 2101 |
2132 |
/** \f[ y_u + y_v + \sum_{B \in \mathcal{O}, uv \in \gamma(B)}z_B \ge
|
| 2102 |
2133 |
w_{uv} \quad \forall uv\in E\f] */
|
| 2103 |
2134 |
/// \f[z_B \ge 0 \quad \forall B \in \mathcal{O}\f]
|
| 2104 |
2135 |
/** \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}}
|
| 2105 |
2136 |
\frac{\vert B \vert - 1}{2}z_B\f] */
|
| 2106 |
2137 |
///
|
| 2107 |
2138 |
/// The algorithm can be executed with the run() function.
|
| 2108 |
2139 |
/// After it the matching (the primal solution) and the dual solution
|
| 2109 |
2140 |
/// can be obtained using the query functions and the
|
| 2110 |
2141 |
/// \ref MaxWeightedPerfectMatching::BlossomIt "BlossomIt" nested class,
|
| 2111 |
2142 |
/// which is able to iterate on the nodes of a blossom.
|
| 2112 |
2143 |
/// If the value type is integer, then the dual solution is multiplied
|
| 2113 |
2144 |
/// by \ref MaxWeightedMatching::dualScale "4".
|
| 2114 |
2145 |
///
|
| 2115 |
2146 |
/// \tparam GR The undirected graph type the algorithm runs on.
|
| 2116 |
2147 |
/// \tparam WM The type edge weight map. The default type is
|
| 2117 |
2148 |
/// \ref concepts::Graph::EdgeMap "GR::EdgeMap<int>".
|
| 2118 |
2149 |
#ifdef DOXYGEN
|
| 2119 |
2150 |
template <typename GR, typename WM>
|
| 2120 |
2151 |
#else
|
| 2121 |
2152 |
template <typename GR,
|
| 2122 |
2153 |
typename WM = typename GR::template EdgeMap<int> >
|
| 2123 |
2154 |
#endif
|
| 2124 |
2155 |
class MaxWeightedPerfectMatching {
|
| 2125 |
2156 |
public:
|
| 2126 |
2157 |
|
| 2127 |
2158 |
/// The graph type of the algorithm
|
| 2128 |
2159 |
typedef GR Graph;
|
| 2129 |
2160 |
/// The type of the edge weight map
|
| 2130 |
2161 |
typedef WM WeightMap;
|
| 2131 |
2162 |
/// The value type of the edge weights
|
| 2132 |
2163 |
typedef typename WeightMap::Value Value;
|
| 2133 |
2164 |
|
| 2134 |
2165 |
/// \brief Scaling factor for dual solution
|
| 2135 |
2166 |
///
|
| 2136 |
2167 |
/// Scaling factor for dual solution, it is equal to 4 or 1
|
| 2137 |
2168 |
/// according to the value type.
|
| 2138 |
2169 |
static const int dualScale =
|
| 2139 |
2170 |
std::numeric_limits<Value>::is_integer ? 4 : 1;
|
| 2140 |
2171 |
|
| 2141 |
2172 |
/// The type of the matching map
|
| 2142 |
2173 |
typedef typename Graph::template NodeMap<typename Graph::Arc>
|
| 2143 |
2174 |
MatchingMap;
|
| 2144 |
2175 |
|
| 2145 |
2176 |
private:
|
| 2146 |
2177 |
|
| 2147 |
2178 |
TEMPLATE_GRAPH_TYPEDEFS(Graph);
|
| 2148 |
2179 |
|
| 2149 |
2180 |
typedef typename Graph::template NodeMap<Value> NodePotential;
|
| 2150 |
2181 |
typedef std::vector<Node> BlossomNodeList;
|
| 2151 |
2182 |
|
| 2152 |
2183 |
struct BlossomVariable {
|
| 2153 |
2184 |
int begin, end;
|
| 2154 |
2185 |
Value value;
|
| 2155 |
2186 |
|
| 2156 |
2187 |
BlossomVariable(int _begin, int _end, Value _value)
|
| 2157 |
2188 |
: begin(_begin), end(_end), value(_value) {}
|
| 2158 |
2189 |
|
| 2159 |
2190 |
};
|
| 2160 |
2191 |
|
| 2161 |
2192 |
typedef std::vector<BlossomVariable> BlossomPotential;
|
| 2162 |
2193 |
|
| 2163 |
2194 |
const Graph& _graph;
|
| 2164 |
2195 |
const WeightMap& _weight;
|
| 2165 |
2196 |
|
| 2166 |
2197 |
MatchingMap* _matching;
|
| 2167 |
2198 |
|
| 2168 |
2199 |
NodePotential* _node_potential;
|
| 2169 |
2200 |
|
| 2170 |
2201 |
BlossomPotential _blossom_potential;
|
| 2171 |
2202 |
BlossomNodeList _blossom_node_list;
|
| 2172 |
2203 |
|
| 2173 |
2204 |
int _node_num;
|
| 2174 |
2205 |
int _blossom_num;
|
| 2175 |
2206 |
|
| 2176 |
2207 |
typedef RangeMap<int> IntIntMap;
|
| 2177 |
2208 |
|
| 2178 |
2209 |
enum Status {
|
| 2179 |
2210 |
EVEN = -1, MATCHED = 0, ODD = 1
|
| 2180 |
2211 |
};
|
| 2181 |
2212 |
|
| 2182 |
2213 |
typedef HeapUnionFind<Value, IntNodeMap> BlossomSet;
|
| 2183 |
2214 |
struct BlossomData {
|
| 2184 |
2215 |
int tree;
|
| 2185 |
2216 |
Status status;
|
| 2186 |
2217 |
Arc pred, next;
|
| 2187 |
2218 |
Value pot, offset;
|
| 2188 |
2219 |
};
|
| 2189 |
2220 |
|
| 2190 |
2221 |
IntNodeMap *_blossom_index;
|
| 2191 |
2222 |
BlossomSet *_blossom_set;
|
| 2192 |
2223 |
RangeMap<BlossomData>* _blossom_data;
|
| 2193 |
2224 |
|
| 2194 |
2225 |
IntNodeMap *_node_index;
|
| 2195 |
2226 |
IntArcMap *_node_heap_index;
|
| 2196 |
2227 |
|
| 2197 |
2228 |
struct NodeData {
|
| 2198 |
2229 |
|
| 2199 |
2230 |
NodeData(IntArcMap& node_heap_index)
|
| 2200 |
2231 |
: heap(node_heap_index) {}
|
| 2201 |
2232 |
|
| 2202 |
2233 |
int blossom;
|
| 2203 |
2234 |
Value pot;
|
| 2204 |
2235 |
BinHeap<Value, IntArcMap> heap;
|
| 2205 |
2236 |
std::map<int, Arc> heap_index;
|
| 2206 |
2237 |
|
| 2207 |
2238 |
int tree;
|
| 2208 |
2239 |
};
|
| 2209 |
2240 |
|
| 2210 |
2241 |
RangeMap<NodeData>* _node_data;
|
| 2211 |
2242 |
|
| 2212 |
2243 |
typedef ExtendFindEnum<IntIntMap> TreeSet;
|
| 2213 |
2244 |
|
| 2214 |
2245 |
IntIntMap *_tree_set_index;
|
| 2215 |
2246 |
TreeSet *_tree_set;
|
| 2216 |
2247 |
|
| 2217 |
2248 |
IntIntMap *_delta2_index;
|
| 2218 |
2249 |
BinHeap<Value, IntIntMap> *_delta2;
|
| 2219 |
2250 |
|
| 2220 |
2251 |
IntEdgeMap *_delta3_index;
|
| 2221 |
2252 |
BinHeap<Value, IntEdgeMap> *_delta3;
|
| 2222 |
2253 |
|
| 2223 |
2254 |
IntIntMap *_delta4_index;
|
| 2224 |
2255 |
BinHeap<Value, IntIntMap> *_delta4;
|
| 2225 |
2256 |
|
| 2226 |
2257 |
Value _delta_sum;
|
| 2227 |
2258 |
int _unmatched;
|
| 2228 |
2259 |
|
| 2229 |
2260 |
typedef MaxWeightedPerfectFractionalMatching<Graph, WeightMap>
|
| 2230 |
2261 |
FractionalMatching;
|
| 2231 |
2262 |
FractionalMatching *_fractional;
|
| 2232 |
2263 |
|
| 2233 |
2264 |
void createStructures() {
|
| 2234 |
2265 |
_node_num = countNodes(_graph);
|
| 2235 |
2266 |
_blossom_num = _node_num * 3 / 2;
|
| 2236 |
2267 |
|
| 2237 |
2268 |
if (!_matching) {
|
| 2238 |
2269 |
_matching = new MatchingMap(_graph);
|
| 2239 |
2270 |
}
|
|
2271 |
|
| 2240 |
2272 |
if (!_node_potential) {
|
| 2241 |
2273 |
_node_potential = new NodePotential(_graph);
|
| 2242 |
2274 |
}
|
|
2275 |
|
| 2243 |
2276 |
if (!_blossom_set) {
|
| 2244 |
2277 |
_blossom_index = new IntNodeMap(_graph);
|
| 2245 |
2278 |
_blossom_set = new BlossomSet(*_blossom_index);
|
| 2246 |
2279 |
_blossom_data = new RangeMap<BlossomData>(_blossom_num);
|
|
2280 |
} else if (_blossom_data->size() != _blossom_num) {
|
|
2281 |
delete _blossom_data;
|
|
2282 |
_blossom_data = new RangeMap<BlossomData>(_blossom_num);
|
| 2247 |
2283 |
}
|
| 2248 |
2284 |
|
| 2249 |
2285 |
if (!_node_index) {
|
| 2250 |
2286 |
_node_index = new IntNodeMap(_graph);
|
| 2251 |
2287 |
_node_heap_index = new IntArcMap(_graph);
|
| 2252 |
2288 |
_node_data = new RangeMap<NodeData>(_node_num,
|
| 2253 |
2289 |
NodeData(*_node_heap_index));
|
|
2290 |
} else if (_node_data->size() != _node_num) {
|
|
2291 |
delete _node_data;
|
|
2292 |
_node_data = new RangeMap<NodeData>(_node_num,
|
|
2293 |
NodeData(*_node_heap_index));
|
| 2254 |
2294 |
}
|
| 2255 |
2295 |
|
| 2256 |
2296 |
if (!_tree_set) {
|
| 2257 |
2297 |
_tree_set_index = new IntIntMap(_blossom_num);
|
| 2258 |
2298 |
_tree_set = new TreeSet(*_tree_set_index);
|
|
2299 |
} else {
|
|
2300 |
_tree_set_index->resize(_blossom_num);
|
| 2259 |
2301 |
}
|
|
2302 |
|
| 2260 |
2303 |
if (!_delta2) {
|
| 2261 |
2304 |
_delta2_index = new IntIntMap(_blossom_num);
|
| 2262 |
2305 |
_delta2 = new BinHeap<Value, IntIntMap>(*_delta2_index);
|
|
2306 |
} else {
|
|
2307 |
_delta2_index->resize(_blossom_num);
|
| 2263 |
2308 |
}
|
|
2309 |
|
| 2264 |
2310 |
if (!_delta3) {
|
| 2265 |
2311 |
_delta3_index = new IntEdgeMap(_graph);
|
| 2266 |
2312 |
_delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index);
|
| 2267 |
2313 |
}
|
|
2314 |
|
| 2268 |
2315 |
if (!_delta4) {
|
| 2269 |
2316 |
_delta4_index = new IntIntMap(_blossom_num);
|
| 2270 |
2317 |
_delta4 = new BinHeap<Value, IntIntMap>(*_delta4_index);
|
|
2318 |
} else {
|
|
2319 |
_delta4_index->resize(_blossom_num);
|
| 2271 |
2320 |
}
|
| 2272 |
2321 |
}
|
| 2273 |
2322 |
|
| 2274 |
2323 |
void destroyStructures() {
|
| 2275 |
2324 |
if (_matching) {
|
| 2276 |
2325 |
delete _matching;
|
| 2277 |
2326 |
}
|
| 2278 |
2327 |
if (_node_potential) {
|
| 2279 |
2328 |
delete _node_potential;
|
| 2280 |
2329 |
}
|
| 2281 |
2330 |
if (_blossom_set) {
|
| 2282 |
2331 |
delete _blossom_index;
|
| 2283 |
2332 |
delete _blossom_set;
|
| 2284 |
2333 |
delete _blossom_data;
|
| 2285 |
2334 |
}
|
| 2286 |
2335 |
|
| 2287 |
2336 |
if (_node_index) {
|
| 2288 |
2337 |
delete _node_index;
|
| 2289 |
2338 |
delete _node_heap_index;
|
| 2290 |
2339 |
delete _node_data;
|
| 2291 |
2340 |
}
|
| 2292 |
2341 |
|
| 2293 |
2342 |
if (_tree_set) {
|
| 2294 |
2343 |
delete _tree_set_index;
|
| 2295 |
2344 |
delete _tree_set;
|
| 2296 |
2345 |
}
|
| 2297 |
2346 |
if (_delta2) {
|
| 2298 |
2347 |
delete _delta2_index;
|
| 2299 |
2348 |
delete _delta2;
|
| 2300 |
2349 |
}
|
| 2301 |
2350 |
if (_delta3) {
|
| 2302 |
2351 |
delete _delta3_index;
|
| 2303 |
2352 |
delete _delta3;
|
| 2304 |
2353 |
}
|
| 2305 |
2354 |
if (_delta4) {
|
| 2306 |
2355 |
delete _delta4_index;
|
| 2307 |
2356 |
delete _delta4;
|
| 2308 |
2357 |
}
|
| 2309 |
2358 |
}
|
| 2310 |
2359 |
|
| 2311 |
2360 |
void matchedToEven(int blossom, int tree) {
|
| 2312 |
2361 |
if (_delta2->state(blossom) == _delta2->IN_HEAP) {
|
| 2313 |
2362 |
_delta2->erase(blossom);
|
| 2314 |
2363 |
}
|
| 2315 |
2364 |
|
| 2316 |
2365 |
if (!_blossom_set->trivial(blossom)) {
|
| 2317 |
2366 |
(*_blossom_data)[blossom].pot -=
|
| 2318 |
2367 |
2 * (_delta_sum - (*_blossom_data)[blossom].offset);
|
| 2319 |
2368 |
}
|
| 2320 |
2369 |
|
| 2321 |
2370 |
for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
|
| 2322 |
2371 |
n != INVALID; ++n) {
|
| 2323 |
2372 |
|
| 2324 |
2373 |
_blossom_set->increase(n, std::numeric_limits<Value>::max());
|
| 2325 |
2374 |
int ni = (*_node_index)[n];
|
| 2326 |
2375 |
|
| 2327 |
2376 |
(*_node_data)[ni].heap.clear();
|
| 2328 |
2377 |
(*_node_data)[ni].heap_index.clear();
|
| 2329 |
2378 |
|
| 2330 |
2379 |
(*_node_data)[ni].pot += _delta_sum - (*_blossom_data)[blossom].offset;
|
| 2331 |
2380 |
|
| 2332 |
2381 |
for (InArcIt e(_graph, n); e != INVALID; ++e) {
|
| 2333 |
2382 |
Node v = _graph.source(e);
|
| 2334 |
2383 |
int vb = _blossom_set->find(v);
|
| 2335 |
2384 |
int vi = (*_node_index)[v];
|
| 2336 |
2385 |
|
| 2337 |
2386 |
Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
|
| 2338 |
2387 |
dualScale * _weight[e];
|
| 2339 |
2388 |
|
| 2340 |
2389 |
if ((*_blossom_data)[vb].status == EVEN) {
|
| 2341 |
2390 |
if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
|
| 2342 |
2391 |
_delta3->push(e, rw / 2);
|
| 2343 |
2392 |
}
|
| 2344 |
2393 |
} else {
|
| 2345 |
2394 |
typename std::map<int, Arc>::iterator it =
|
| 2346 |
2395 |
(*_node_data)[vi].heap_index.find(tree);
|
| 2347 |
2396 |
|
| 2348 |
2397 |
if (it != (*_node_data)[vi].heap_index.end()) {
|
| 2349 |
2398 |
if ((*_node_data)[vi].heap[it->second] > rw) {
|
| 2350 |
2399 |
(*_node_data)[vi].heap.replace(it->second, e);
|
| 2351 |
2400 |
(*_node_data)[vi].heap.decrease(e, rw);
|
| 2352 |
2401 |
it->second = e;
|
| 2353 |
2402 |
}
|
| 2354 |
2403 |
} else {
|
| 2355 |
2404 |
(*_node_data)[vi].heap.push(e, rw);
|
| 2356 |
2405 |
(*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
|
| 2357 |
2406 |
}
|
| 2358 |
2407 |
|
| 2359 |
2408 |
if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
|
| 2360 |
2409 |
_blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
|
| 2361 |
2410 |
|
| 2362 |
2411 |
if ((*_blossom_data)[vb].status == MATCHED) {
|
| 2363 |
2412 |
if (_delta2->state(vb) != _delta2->IN_HEAP) {
|
| 2364 |
2413 |
_delta2->push(vb, _blossom_set->classPrio(vb) -
|
| 2365 |
2414 |
(*_blossom_data)[vb].offset);
|
| 2366 |
2415 |
} else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
|
| 2367 |
2416 |
(*_blossom_data)[vb].offset){
|
| 2368 |
2417 |
_delta2->decrease(vb, _blossom_set->classPrio(vb) -
|
| 2369 |
2418 |
(*_blossom_data)[vb].offset);
|
| 2370 |
2419 |
}
|
| 2371 |
2420 |
}
|
| 2372 |
2421 |
}
|
| 2373 |
2422 |
}
|
| 2374 |
2423 |
}
|
| 2375 |
2424 |
}
|
| 2376 |
2425 |
(*_blossom_data)[blossom].offset = 0;
|
| 2377 |
2426 |
}
|
| 2378 |
2427 |
|
| 2379 |
2428 |
void matchedToOdd(int blossom) {
|
| 2380 |
2429 |
if (_delta2->state(blossom) == _delta2->IN_HEAP) {
|
| 2381 |
2430 |
_delta2->erase(blossom);
|
| 2382 |
2431 |
}
|
| 2383 |
2432 |
(*_blossom_data)[blossom].offset += _delta_sum;
|
| 2384 |
2433 |
if (!_blossom_set->trivial(blossom)) {
|
| 2385 |
2434 |
_delta4->push(blossom, (*_blossom_data)[blossom].pot / 2 +
|
| 2386 |
2435 |
(*_blossom_data)[blossom].offset);
|
| 2387 |
2436 |
}
|
| 2388 |
2437 |
}
|
| 2389 |
2438 |
|
| 2390 |
2439 |
void evenToMatched(int blossom, int tree) {
|
| 2391 |
2440 |
if (!_blossom_set->trivial(blossom)) {
|
| 2392 |
2441 |
(*_blossom_data)[blossom].pot += 2 * _delta_sum;
|
| 2393 |
2442 |
}
|
| 2394 |
2443 |
|
| 2395 |
2444 |
for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
|
| 2396 |
2445 |
n != INVALID; ++n) {
|
| 2397 |
2446 |
int ni = (*_node_index)[n];
|
| 2398 |
2447 |
(*_node_data)[ni].pot -= _delta_sum;
|
| 2399 |
2448 |
|
| 2400 |
2449 |
for (InArcIt e(_graph, n); e != INVALID; ++e) {
|
| 2401 |
2450 |
Node v = _graph.source(e);
|
| 2402 |
2451 |
int vb = _blossom_set->find(v);
|
| 2403 |
2452 |
int vi = (*_node_index)[v];
|
| 2404 |
2453 |
|
| 2405 |
2454 |
Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
|
| 2406 |
2455 |
dualScale * _weight[e];
|
| 2407 |
2456 |
|
| 2408 |
2457 |
if (vb == blossom) {
|
| 2409 |
2458 |
if (_delta3->state(e) == _delta3->IN_HEAP) {
|
| 2410 |
2459 |
_delta3->erase(e);
|
| 2411 |
2460 |
}
|
| 2412 |
2461 |
} else if ((*_blossom_data)[vb].status == EVEN) {
|
| 2413 |
2462 |
|
| 2414 |
2463 |
if (_delta3->state(e) == _delta3->IN_HEAP) {
|
| 2415 |
2464 |
_delta3->erase(e);
|
| 2416 |
2465 |
}
|
| 2417 |
2466 |
|
| 2418 |
2467 |
int vt = _tree_set->find(vb);
|
| 2419 |
2468 |
|
| 2420 |
2469 |
if (vt != tree) {
|
| 2421 |
2470 |
|
| 2422 |
2471 |
Arc r = _graph.oppositeArc(e);
|
| 2423 |
2472 |
|
| 2424 |
2473 |
typename std::map<int, Arc>::iterator it =
|
| 2425 |
2474 |
(*_node_data)[ni].heap_index.find(vt);
|
| 2426 |
2475 |
|
| 2427 |
2476 |
if (it != (*_node_data)[ni].heap_index.end()) {
|
| 2428 |
2477 |
if ((*_node_data)[ni].heap[it->second] > rw) {
|
| 2429 |
2478 |
(*_node_data)[ni].heap.replace(it->second, r);
|
| 2430 |
2479 |
(*_node_data)[ni].heap.decrease(r, rw);
|
| 2431 |
2480 |
it->second = r;
|
| 2432 |
2481 |
}
|
| 2433 |
2482 |
} else {
|
| 2434 |
2483 |
(*_node_data)[ni].heap.push(r, rw);
|
| 2435 |
2484 |
(*_node_data)[ni].heap_index.insert(std::make_pair(vt, r));
|
| 2436 |
2485 |
}
|
| 2437 |
2486 |
|
| 2438 |
2487 |
if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) {
|
| 2439 |
2488 |
_blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
|
| 2440 |
2489 |
|
| 2441 |
2490 |
if (_delta2->state(blossom) != _delta2->IN_HEAP) {
|
| 2442 |
2491 |
_delta2->push(blossom, _blossom_set->classPrio(blossom) -
|
| 2443 |
2492 |
(*_blossom_data)[blossom].offset);
|
| 2444 |
2493 |
} else if ((*_delta2)[blossom] >
|
| 2445 |
2494 |
_blossom_set->classPrio(blossom) -
|
| 2446 |
2495 |
(*_blossom_data)[blossom].offset){
|
| 2447 |
2496 |
_delta2->decrease(blossom, _blossom_set->classPrio(blossom) -
|
| 2448 |
2497 |
(*_blossom_data)[blossom].offset);
|
| 2449 |
2498 |
}
|
| 2450 |
2499 |
}
|
| 2451 |
2500 |
}
|
| 2452 |
2501 |
} else {
|
| 2453 |
2502 |
|
| 2454 |
2503 |
typename std::map<int, Arc>::iterator it =
|
| 2455 |
2504 |
(*_node_data)[vi].heap_index.find(tree);
|
| 2456 |
2505 |
|
| 2457 |
2506 |
if (it != (*_node_data)[vi].heap_index.end()) {
|
| 2458 |
2507 |
(*_node_data)[vi].heap.erase(it->second);
|
| 2459 |
2508 |
(*_node_data)[vi].heap_index.erase(it);
|
| 2460 |
2509 |
if ((*_node_data)[vi].heap.empty()) {
|
| 2461 |
2510 |
_blossom_set->increase(v, std::numeric_limits<Value>::max());
|
| 2462 |
2511 |
} else if ((*_blossom_set)[v] < (*_node_data)[vi].heap.prio()) {
|
| 2463 |
2512 |
_blossom_set->increase(v, (*_node_data)[vi].heap.prio());
|
| 2464 |
2513 |
}
|
| 2465 |
2514 |
|
| 2466 |
2515 |
if ((*_blossom_data)[vb].status == MATCHED) {
|
| 2467 |
2516 |
if (_blossom_set->classPrio(vb) ==
|
| 2468 |
2517 |
std::numeric_limits<Value>::max()) {
|
| 2469 |
2518 |
_delta2->erase(vb);
|
| 2470 |
2519 |
} else if ((*_delta2)[vb] < _blossom_set->classPrio(vb) -
|
| 2471 |
2520 |
(*_blossom_data)[vb].offset) {
|
| 2472 |
2521 |
_delta2->increase(vb, _blossom_set->classPrio(vb) -
|
| 2473 |
2522 |
(*_blossom_data)[vb].offset);
|
| 2474 |
2523 |
}
|
| 2475 |
2524 |
}
|
| 2476 |
2525 |
}
|
| 2477 |
2526 |
}
|
| 2478 |
2527 |
}
|
| 2479 |
2528 |
}
|
| 2480 |
2529 |
}
|
| 2481 |
2530 |
|
| 2482 |
2531 |
void oddToMatched(int blossom) {
|
| 2483 |
2532 |
(*_blossom_data)[blossom].offset -= _delta_sum;
|
| 2484 |
2533 |
|
| 2485 |
2534 |
if (_blossom_set->classPrio(blossom) !=
|
| 2486 |
2535 |
std::numeric_limits<Value>::max()) {
|
| 2487 |
2536 |
_delta2->push(blossom, _blossom_set->classPrio(blossom) -
|
| 2488 |
2537 |
(*_blossom_data)[blossom].offset);
|
| 2489 |
2538 |
}
|
| 2490 |
2539 |
|
| 2491 |
2540 |
if (!_blossom_set->trivial(blossom)) {
|
| 2492 |
2541 |
_delta4->erase(blossom);
|
| 2493 |
2542 |
}
|
| 2494 |
2543 |
}
|
| 2495 |
2544 |
|
| 2496 |
2545 |
void oddToEven(int blossom, int tree) {
|
| 2497 |
2546 |
if (!_blossom_set->trivial(blossom)) {
|
| 2498 |
2547 |
_delta4->erase(blossom);
|
| 2499 |
2548 |
(*_blossom_data)[blossom].pot -=
|
| 2500 |
2549 |
2 * (2 * _delta_sum - (*_blossom_data)[blossom].offset);
|
| 2501 |
2550 |
}
|
| 2502 |
2551 |
|
| 2503 |
2552 |
for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
|
| 2504 |
2553 |
n != INVALID; ++n) {
|
| 2505 |
2554 |
int ni = (*_node_index)[n];
|
| 2506 |
2555 |
|
| 2507 |
2556 |
_blossom_set->increase(n, std::numeric_limits<Value>::max());
|
| 2508 |
2557 |
|
| 2509 |
2558 |
(*_node_data)[ni].heap.clear();
|
| 2510 |
2559 |
(*_node_data)[ni].heap_index.clear();
|
| 2511 |
2560 |
(*_node_data)[ni].pot +=
|
| 2512 |
2561 |
2 * _delta_sum - (*_blossom_data)[blossom].offset;
|
| 2513 |
2562 |
|
| 2514 |
2563 |
for (InArcIt e(_graph, n); e != INVALID; ++e) {
|
| 2515 |
2564 |
Node v = _graph.source(e);
|
| 2516 |
2565 |
int vb = _blossom_set->find(v);
|
| 2517 |
2566 |
int vi = (*_node_index)[v];
|
| 2518 |
2567 |
|
| 2519 |
2568 |
Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
|
| 2520 |
2569 |
dualScale * _weight[e];
|
| 2521 |
2570 |
|
| 2522 |
2571 |
if ((*_blossom_data)[vb].status == EVEN) {
|
| 2523 |
2572 |
if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
|
| 2524 |
2573 |
_delta3->push(e, rw / 2);
|
| 2525 |
2574 |
}
|
| 2526 |
2575 |
} else {
|
| 2527 |
2576 |
|
| 2528 |
2577 |
typename std::map<int, Arc>::iterator it =
|
| 2529 |
2578 |
(*_node_data)[vi].heap_index.find(tree);
|
| 2530 |
2579 |
|
| 2531 |
2580 |
if (it != (*_node_data)[vi].heap_index.end()) {
|
| 2532 |
2581 |
if ((*_node_data)[vi].heap[it->second] > rw) {
|
| 2533 |
2582 |
(*_node_data)[vi].heap.replace(it->second, e);
|
| 2534 |
2583 |
(*_node_data)[vi].heap.decrease(e, rw);
|
| 2535 |
2584 |
it->second = e;
|
| 2536 |
2585 |
}
|
| 2537 |
2586 |
} else {
|
| 2538 |
2587 |
(*_node_data)[vi].heap.push(e, rw);
|
| 2539 |
2588 |
(*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
|
| 2540 |
2589 |
}
|
| 2541 |
2590 |
|
| 2542 |
2591 |
if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
|
| 2543 |
2592 |
_blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
|
| 2544 |
2593 |
|
| 2545 |
2594 |
if ((*_blossom_data)[vb].status == MATCHED) {
|
| 2546 |
2595 |
if (_delta2->state(vb) != _delta2->IN_HEAP) {
|
| 2547 |
2596 |
_delta2->push(vb, _blossom_set->classPrio(vb) -
|
| 2548 |
2597 |
(*_blossom_data)[vb].offset);
|
| 2549 |
2598 |
} else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
|
| 2550 |
2599 |
(*_blossom_data)[vb].offset) {
|
| 2551 |
2600 |
_delta2->decrease(vb, _blossom_set->classPrio(vb) -
|
| 2552 |
2601 |
(*_blossom_data)[vb].offset);
|
| 2553 |
2602 |
}
|
| 2554 |
2603 |
}
|
| 2555 |
2604 |
}
|
| 2556 |
2605 |
}
|
| 2557 |
2606 |
}
|
| 2558 |
2607 |
}
|
| 2559 |
2608 |
(*_blossom_data)[blossom].offset = 0;
|
| 2560 |
2609 |
}
|
| 2561 |
2610 |
|
| 2562 |
2611 |
void alternatePath(int even, int tree) {
|
| 2563 |
2612 |
int odd;
|
| 2564 |
2613 |
|
| 2565 |
2614 |
evenToMatched(even, tree);
|
| 2566 |
2615 |
(*_blossom_data)[even].status = MATCHED;
|
| 2567 |
2616 |
|
| 2568 |
2617 |
while ((*_blossom_data)[even].pred != INVALID) {
|
| 2569 |
2618 |
odd = _blossom_set->find(_graph.target((*_blossom_data)[even].pred));
|
| 2570 |
2619 |
(*_blossom_data)[odd].status = MATCHED;
|
| 2571 |
2620 |
oddToMatched(odd);
|
| 2572 |
2621 |
(*_blossom_data)[odd].next = (*_blossom_data)[odd].pred;
|
| 2573 |
2622 |
|
| 2574 |
2623 |
even = _blossom_set->find(_graph.target((*_blossom_data)[odd].pred));
|
| 2575 |
2624 |
(*_blossom_data)[even].status = MATCHED;
|
| 2576 |
2625 |
evenToMatched(even, tree);
|
| 2577 |
2626 |
(*_blossom_data)[even].next =
|
| 2578 |
2627 |
_graph.oppositeArc((*_blossom_data)[odd].pred);
|
| 2579 |
2628 |
}
|
| 2580 |
2629 |
|
| 2581 |
2630 |
}
|
| 2582 |
2631 |
|
| 2583 |
2632 |
void destroyTree(int tree) {
|
| 2584 |
2633 |
for (TreeSet::ItemIt b(*_tree_set, tree); b != INVALID; ++b) {
|
| 2585 |
2634 |
if ((*_blossom_data)[b].status == EVEN) {
|
| 2586 |
2635 |
(*_blossom_data)[b].status = MATCHED;
|
| 2587 |
2636 |
evenToMatched(b, tree);
|
| 2588 |
2637 |
} else if ((*_blossom_data)[b].status == ODD) {
|
| 2589 |
2638 |
(*_blossom_data)[b].status = MATCHED;
|
| 2590 |
2639 |
oddToMatched(b);
|
| 2591 |
2640 |
}
|
| 2592 |
2641 |
}
|
| 2593 |
2642 |
_tree_set->eraseClass(tree);
|
| 2594 |
2643 |
}
|
| 2595 |
2644 |
|
| 2596 |
2645 |
void augmentOnEdge(const Edge& edge) {
|
| 2597 |
2646 |
|
| 2598 |
2647 |
int left = _blossom_set->find(_graph.u(edge));
|
| 2599 |
2648 |
int right = _blossom_set->find(_graph.v(edge));
|
| 2600 |
2649 |
|
| 2601 |
2650 |
int left_tree = _tree_set->find(left);
|
| 2602 |
2651 |
alternatePath(left, left_tree);
|
| 2603 |
2652 |
destroyTree(left_tree);
|
| 2604 |
2653 |
|
| 2605 |
2654 |
int right_tree = _tree_set->find(right);
|
| 2606 |
2655 |
alternatePath(right, right_tree);
|
| 2607 |
2656 |
destroyTree(right_tree);
|
| 2608 |
2657 |
|
| 2609 |
2658 |
(*_blossom_data)[left].next = _graph.direct(edge, true);
|
| 2610 |
2659 |
(*_blossom_data)[right].next = _graph.direct(edge, false);
|
| 2611 |
2660 |
}
|
| 2612 |
2661 |
|
| 2613 |
2662 |
void extendOnArc(const Arc& arc) {
|
| 2614 |
2663 |
int base = _blossom_set->find(_graph.target(arc));
|
| 2615 |
2664 |
int tree = _tree_set->find(base);
|
| 2616 |
2665 |
|
| 2617 |
2666 |
int odd = _blossom_set->find(_graph.source(arc));
|
| 2618 |
2667 |
_tree_set->insert(odd, tree);
|
| 2619 |
2668 |
(*_blossom_data)[odd].status = ODD;
|
| 2620 |
2669 |
matchedToOdd(odd);
|
| 2621 |
2670 |
(*_blossom_data)[odd].pred = arc;
|
| 2622 |
2671 |
|
| 2623 |
2672 |
int even = _blossom_set->find(_graph.target((*_blossom_data)[odd].next));
|
| 2624 |
2673 |
(*_blossom_data)[even].pred = (*_blossom_data)[even].next;
|
| 2625 |
2674 |
_tree_set->insert(even, tree);
|
| 2626 |
2675 |
(*_blossom_data)[even].status = EVEN;
|
| 2627 |
2676 |
matchedToEven(even, tree);
|
| 2628 |
2677 |
}
|
| 2629 |
2678 |
|
| 2630 |
2679 |
void shrinkOnEdge(const Edge& edge, int tree) {
|
| 2631 |
2680 |
int nca = -1;
|
| 2632 |
2681 |
std::vector<int> left_path, right_path;
|
| 2633 |
2682 |
|
| 2634 |
2683 |
{
|
| 2635 |
2684 |
std::set<int> left_set, right_set;
|
| 2636 |
2685 |
int left = _blossom_set->find(_graph.u(edge));
|
| 2637 |
2686 |
left_path.push_back(left);
|
| 2638 |
2687 |
left_set.insert(left);
|
| 2639 |
2688 |
|
| 2640 |
2689 |
int right = _blossom_set->find(_graph.v(edge));
|
| 2641 |
2690 |
right_path.push_back(right);
|
| 2642 |
2691 |
right_set.insert(right);
|
| 2643 |
2692 |
|
| 2644 |
2693 |
while (true) {
|
| 2645 |
2694 |
|
| 2646 |
2695 |
if ((*_blossom_data)[left].pred == INVALID) break;
|
| 2647 |
2696 |
|
| 2648 |
2697 |
left =
|
| 2649 |
2698 |
_blossom_set->find(_graph.target((*_blossom_data)[left].pred));
|
| 2650 |
2699 |
left_path.push_back(left);
|
| 2651 |
2700 |
left =
|
| 2652 |
2701 |
_blossom_set->find(_graph.target((*_blossom_data)[left].pred));
|
| 2653 |
2702 |
left_path.push_back(left);
|
| 2654 |
2703 |
|
| 2655 |
2704 |
left_set.insert(left);
|
| 2656 |
2705 |
|
| 2657 |
2706 |
if (right_set.find(left) != right_set.end()) {
|
| 2658 |
2707 |
nca = left;
|
| 2659 |
2708 |
break;
|
| 2660 |
2709 |
}
|
| 2661 |
2710 |
|
| 2662 |
2711 |
if ((*_blossom_data)[right].pred == INVALID) break;
|
| 2663 |
2712 |
|
| 2664 |
2713 |
right =
|
| 2665 |
2714 |
_blossom_set->find(_graph.target((*_blossom_data)[right].pred));
|
| 2666 |
2715 |
right_path.push_back(right);
|
| 2667 |
2716 |
right =
|
| 2668 |
2717 |
_blossom_set->find(_graph.target((*_blossom_data)[right].pred));
|
| 2669 |
2718 |
right_path.push_back(right);
|
| 2670 |
2719 |
|
| 2671 |
2720 |
right_set.insert(right);
|
| 2672 |
2721 |
|
| 2673 |
2722 |
if (left_set.find(right) != left_set.end()) {
|
| 2674 |
2723 |
nca = right;
|
| 2675 |
2724 |
break;
|
| 2676 |
2725 |
}
|
| 2677 |
2726 |
|
| 2678 |
2727 |
}
|
| 2679 |
2728 |
|
| 2680 |
2729 |
if (nca == -1) {
|
| 2681 |
2730 |
if ((*_blossom_data)[left].pred == INVALID) {
|
| 2682 |
2731 |
nca = right;
|
| 2683 |
2732 |
while (left_set.find(nca) == left_set.end()) {
|
| 2684 |
2733 |
nca =
|
| 2685 |
2734 |
_blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
|
| 2686 |
2735 |
right_path.push_back(nca);
|
| 2687 |
2736 |
nca =
|
| 2688 |
2737 |
_blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
|
| 2689 |
2738 |
right_path.push_back(nca);
|
| 2690 |
2739 |
}
|
| 2691 |
2740 |
} else {
|
| 2692 |
2741 |
nca = left;
|
| 2693 |
2742 |
while (right_set.find(nca) == right_set.end()) {
|
| 2694 |
2743 |
nca =
|
| 2695 |
2744 |
_blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
|
| 2696 |
2745 |
left_path.push_back(nca);
|
| 2697 |
2746 |
nca =
|
| 2698 |
2747 |
_blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
|
| 2699 |
2748 |
left_path.push_back(nca);
|
| 2700 |
2749 |
}
|
| 2701 |
2750 |
}
|
| 2702 |
2751 |
}
|
| 2703 |
2752 |
}
|
| 2704 |
2753 |
|
| 2705 |
2754 |
std::vector<int> subblossoms;
|
| 2706 |
2755 |
Arc prev;
|
| 2707 |
2756 |
|
| 2708 |
2757 |
prev = _graph.direct(edge, true);
|
| 2709 |
2758 |
for (int i = 0; left_path[i] != nca; i += 2) {
|
| 2710 |
2759 |
subblossoms.push_back(left_path[i]);
|
| 2711 |
2760 |
(*_blossom_data)[left_path[i]].next = prev;
|
| 2712 |
2761 |
_tree_set->erase(left_path[i]);
|
| 2713 |
2762 |
|
| 2714 |
2763 |
subblossoms.push_back(left_path[i + 1]);
|
| 2715 |
2764 |
(*_blossom_data)[left_path[i + 1]].status = EVEN;
|
| 2716 |
2765 |
oddToEven(left_path[i + 1], tree);
|
| 2717 |
2766 |
_tree_set->erase(left_path[i + 1]);
|
| 2718 |
2767 |
prev = _graph.oppositeArc((*_blossom_data)[left_path[i + 1]].pred);
|
| 2719 |
2768 |
}
|
| 2720 |
2769 |
|
| 2721 |
2770 |
int k = 0;
|
| 2722 |
2771 |
while (right_path[k] != nca) ++k;
|
| 2723 |
2772 |
|
| 2724 |
2773 |
subblossoms.push_back(nca);
|
| 2725 |
2774 |
(*_blossom_data)[nca].next = prev;
|
| 2726 |
2775 |
|
| 2727 |
2776 |
for (int i = k - 2; i >= 0; i -= 2) {
|
| 2728 |
2777 |
subblossoms.push_back(right_path[i + 1]);
|
| 2729 |
2778 |
(*_blossom_data)[right_path[i + 1]].status = EVEN;
|
| 2730 |
2779 |
oddToEven(right_path[i + 1], tree);
|
| 2731 |
2780 |
_tree_set->erase(right_path[i + 1]);
|
| 2732 |
2781 |
|
| 2733 |
2782 |
(*_blossom_data)[right_path[i + 1]].next =
|
| 2734 |
2783 |
(*_blossom_data)[right_path[i + 1]].pred;
|
| 2735 |
2784 |
|
| 2736 |
2785 |
subblossoms.push_back(right_path[i]);
|
| 2737 |
2786 |
_tree_set->erase(right_path[i]);
|
| 2738 |
2787 |
}
|
| 2739 |
2788 |
|
| 2740 |
2789 |
int surface =
|
| 2741 |
2790 |
_blossom_set->join(subblossoms.begin(), subblossoms.end());
|
| 2742 |
2791 |
|
| 2743 |
2792 |
for (int i = 0; i < int(subblossoms.size()); ++i) {
|
| 2744 |
2793 |
if (!_blossom_set->trivial(subblossoms[i])) {
|
| 2745 |
2794 |
(*_blossom_data)[subblossoms[i]].pot += 2 * _delta_sum;
|
| 2746 |
2795 |
}
|
| 2747 |
2796 |
(*_blossom_data)[subblossoms[i]].status = MATCHED;
|
| 2748 |
2797 |
}
|
| 2749 |
2798 |
|
| 2750 |
2799 |
(*_blossom_data)[surface].pot = -2 * _delta_sum;
|
| 2751 |
2800 |
(*_blossom_data)[surface].offset = 0;
|
| 2752 |
2801 |
(*_blossom_data)[surface].status = EVEN;
|
| 2753 |
2802 |
(*_blossom_data)[surface].pred = (*_blossom_data)[nca].pred;
|
| 2754 |
2803 |
(*_blossom_data)[surface].next = (*_blossom_data)[nca].pred;
|
| 2755 |
2804 |
|
| 2756 |
2805 |
_tree_set->insert(surface, tree);
|
| 2757 |
2806 |
_tree_set->erase(nca);
|
| 2758 |
2807 |
}
|
| 2759 |
2808 |
|
| 2760 |
2809 |
void splitBlossom(int blossom) {
|
| 2761 |
2810 |
Arc next = (*_blossom_data)[blossom].next;
|
| 2762 |
2811 |
Arc pred = (*_blossom_data)[blossom].pred;
|
| 2763 |
2812 |
|
| 2764 |
2813 |
int tree = _tree_set->find(blossom);
|
| 2765 |
2814 |
|
| 2766 |
2815 |
(*_blossom_data)[blossom].status = MATCHED;
|
| 2767 |
2816 |
oddToMatched(blossom);
|
| 2768 |
2817 |
if (_delta2->state(blossom) == _delta2->IN_HEAP) {
|
| 2769 |
2818 |
_delta2->erase(blossom);
|
| 2770 |
2819 |
}
|
| 2771 |
2820 |
|
| 2772 |
2821 |
std::vector<int> subblossoms;
|
| 2773 |
2822 |
_blossom_set->split(blossom, std::back_inserter(subblossoms));
|
| 2774 |
2823 |
|
| 2775 |
2824 |
Value offset = (*_blossom_data)[blossom].offset;
|
| 2776 |
2825 |
int b = _blossom_set->find(_graph.source(pred));
|
| 2777 |
2826 |
int d = _blossom_set->find(_graph.source(next));
|
| 2778 |
2827 |
|
| 2779 |
2828 |
int ib = -1, id = -1;
|
| 2780 |
2829 |
for (int i = 0; i < int(subblossoms.size()); ++i) {
|
| 2781 |
2830 |
if (subblossoms[i] == b) ib = i;
|
| 2782 |
2831 |
if (subblossoms[i] == d) id = i;
|
| 2783 |
2832 |
|
| 2784 |
2833 |
(*_blossom_data)[subblossoms[i]].offset = offset;
|
| 2785 |
2834 |
if (!_blossom_set->trivial(subblossoms[i])) {
|
| 2786 |
2835 |
(*_blossom_data)[subblossoms[i]].pot -= 2 * offset;
|
| 2787 |
2836 |
}
|
| 2788 |
2837 |
if (_blossom_set->classPrio(subblossoms[i]) !=
|
| 2789 |
2838 |
std::numeric_limits<Value>::max()) {
|
| 2790 |
2839 |
_delta2->push(subblossoms[i],
|
| 2791 |
2840 |
_blossom_set->classPrio(subblossoms[i]) -
|
| 2792 |
2841 |
(*_blossom_data)[subblossoms[i]].offset);
|
| 2793 |
2842 |
}
|
| 2794 |
2843 |
}
|
| 2795 |
2844 |
|
| 2796 |
2845 |
if (id > ib ? ((id - ib) % 2 == 0) : ((ib - id) % 2 == 1)) {
|
| 2797 |
2846 |
for (int i = (id + 1) % subblossoms.size();
|
| 2798 |
2847 |
i != ib; i = (i + 2) % subblossoms.size()) {
|
| 2799 |
2848 |
int sb = subblossoms[i];
|
| 2800 |
2849 |
int tb = subblossoms[(i + 1) % subblossoms.size()];
|
| 2801 |
2850 |
(*_blossom_data)[sb].next =
|
| 2802 |
2851 |
_graph.oppositeArc((*_blossom_data)[tb].next);
|
| 2803 |
2852 |
}
|
| 2804 |
2853 |
|
| 2805 |
2854 |
for (int i = ib; i != id; i = (i + 2) % subblossoms.size()) {
|
| 2806 |
2855 |
int sb = subblossoms[i];
|
| 2807 |
2856 |
int tb = subblossoms[(i + 1) % subblossoms.size()];
|
| 2808 |
2857 |
int ub = subblossoms[(i + 2) % subblossoms.size()];
|
| 2809 |
2858 |
|
| 2810 |
2859 |
(*_blossom_data)[sb].status = ODD;
|
| 2811 |
2860 |
matchedToOdd(sb);
|
| 2812 |
2861 |
_tree_set->insert(sb, tree);
|
| 2813 |
2862 |
(*_blossom_data)[sb].pred = pred;
|
| 2814 |
2863 |
(*_blossom_data)[sb].next =
|
| 2815 |
2864 |
_graph.oppositeArc((*_blossom_data)[tb].next);
|
| 2816 |
2865 |
|
| 2817 |
2866 |
pred = (*_blossom_data)[ub].next;
|
| 2818 |
2867 |
|
| 2819 |
2868 |
(*_blossom_data)[tb].status = EVEN;
|
| 2820 |
2869 |
matchedToEven(tb, tree);
|
| 2821 |
2870 |
_tree_set->insert(tb, tree);
|
| 2822 |
2871 |
(*_blossom_data)[tb].pred = (*_blossom_data)[tb].next;
|
| 2823 |
2872 |
}
|
| 2824 |
2873 |
|
| 2825 |
2874 |
(*_blossom_data)[subblossoms[id]].status = ODD;
|
| 2826 |
2875 |
matchedToOdd(subblossoms[id]);
|
| 2827 |
2876 |
_tree_set->insert(subblossoms[id], tree);
|
| 2828 |
2877 |
(*_blossom_data)[subblossoms[id]].next = next;
|
| 2829 |
2878 |
(*_blossom_data)[subblossoms[id]].pred = pred;
|
| 2830 |
2879 |
|
| 2831 |
2880 |
} else {
|
| 2832 |
2881 |
|
| 2833 |
2882 |
for (int i = (ib + 1) % subblossoms.size();
|
| 2834 |
2883 |
i != id; i = (i + 2) % subblossoms.size()) {
|
| 2835 |
2884 |
int sb = subblossoms[i];
|
| 2836 |
2885 |
int tb = subblossoms[(i + 1) % subblossoms.size()];
|
| 2837 |
2886 |
(*_blossom_data)[sb].next =
|
| 2838 |
2887 |
_graph.oppositeArc((*_blossom_data)[tb].next);
|
| 2839 |
2888 |
}
|
| 2840 |
2889 |
|
| 2841 |
2890 |
for (int i = id; i != ib; i = (i + 2) % subblossoms.size()) {
|
| 2842 |
2891 |
int sb = subblossoms[i];
|
| 2843 |
2892 |
int tb = subblossoms[(i + 1) % subblossoms.size()];
|
| 2844 |
2893 |
int ub = subblossoms[(i + 2) % subblossoms.size()];
|
| 2845 |
2894 |
|
| 2846 |
2895 |
(*_blossom_data)[sb].status = ODD;
|
| 2847 |
2896 |
matchedToOdd(sb);
|
| 2848 |
2897 |
_tree_set->insert(sb, tree);
|
| 2849 |
2898 |
(*_blossom_data)[sb].next = next;
|
| 2850 |
2899 |
(*_blossom_data)[sb].pred =
|
| 2851 |
2900 |
_graph.oppositeArc((*_blossom_data)[tb].next);
|
| 2852 |
2901 |
|
| 2853 |
2902 |
(*_blossom_data)[tb].status = EVEN;
|
| 2854 |
2903 |
matchedToEven(tb, tree);
|
| 2855 |
2904 |
_tree_set->insert(tb, tree);
|
| 2856 |
2905 |
(*_blossom_data)[tb].pred =
|
| 2857 |
2906 |
(*_blossom_data)[tb].next =
|
| 2858 |
2907 |
_graph.oppositeArc((*_blossom_data)[ub].next);
|
| 2859 |
2908 |
next = (*_blossom_data)[ub].next;
|
| 2860 |
2909 |
}
|
| 2861 |
2910 |
|
| 2862 |
2911 |
(*_blossom_data)[subblossoms[ib]].status = ODD;
|
| 2863 |
2912 |
matchedToOdd(subblossoms[ib]);
|
| 2864 |
2913 |
_tree_set->insert(subblossoms[ib], tree);
|
| 2865 |
2914 |
(*_blossom_data)[subblossoms[ib]].next = next;
|
| 2866 |
2915 |
(*_blossom_data)[subblossoms[ib]].pred = pred;
|
| 2867 |
2916 |
}
|
| 2868 |
2917 |
_tree_set->erase(blossom);
|
| 2869 |
2918 |
}
|
| 2870 |
2919 |
|
| 2871 |
2920 |
void extractBlossom(int blossom, const Node& base, const Arc& matching) {
|
| 2872 |
2921 |
if (_blossom_set->trivial(blossom)) {
|
| 2873 |
2922 |
int bi = (*_node_index)[base];
|
| 2874 |
2923 |
Value pot = (*_node_data)[bi].pot;
|
| 2875 |
2924 |
|
| 2876 |
2925 |
(*_matching)[base] = matching;
|
| 2877 |
2926 |
_blossom_node_list.push_back(base);
|
| 2878 |
2927 |
(*_node_potential)[base] = pot;
|
| 2879 |
2928 |
} else {
|
| 2880 |
2929 |
|
| 2881 |
2930 |
Value pot = (*_blossom_data)[blossom].pot;
|
| 2882 |
2931 |
int bn = _blossom_node_list.size();
|
| 2883 |
2932 |
|
| 2884 |
2933 |
std::vector<int> subblossoms;
|
| 2885 |
2934 |
_blossom_set->split(blossom, std::back_inserter(subblossoms));
|
| 2886 |
2935 |
int b = _blossom_set->find(base);
|
| 2887 |
2936 |
int ib = -1;
|
| 2888 |
2937 |
for (int i = 0; i < int(subblossoms.size()); ++i) {
|
| 2889 |
2938 |
if (subblossoms[i] == b) { ib = i; break; }
|
| 2890 |
2939 |
}
|
| 2891 |
2940 |
|
| 2892 |
2941 |
for (int i = 1; i < int(subblossoms.size()); i += 2) {
|
| 2893 |
2942 |
int sb = subblossoms[(ib + i) % subblossoms.size()];
|
| 2894 |
2943 |
int tb = subblossoms[(ib + i + 1) % subblossoms.size()];
|
| 2895 |
2944 |
|
| 2896 |
2945 |
Arc m = (*_blossom_data)[tb].next;
|
| 2897 |
2946 |
extractBlossom(sb, _graph.target(m), _graph.oppositeArc(m));
|
| 2898 |
2947 |
extractBlossom(tb, _graph.source(m), m);
|
| 2899 |
2948 |
}
|
| 2900 |
2949 |
extractBlossom(subblossoms[ib], base, matching);
|
| 2901 |
2950 |
|
| 2902 |
2951 |
int en = _blossom_node_list.size();
|
| 2903 |
2952 |
|
| 2904 |
2953 |
_blossom_potential.push_back(BlossomVariable(bn, en, pot));
|
| 2905 |
2954 |
}
|
| 2906 |
2955 |
}
|
| 2907 |
2956 |
|
| 2908 |
2957 |
void extractMatching() {
|
| 2909 |
2958 |
std::vector<int> blossoms;
|
| 2910 |
2959 |
for (typename BlossomSet::ClassIt c(*_blossom_set); c != INVALID; ++c) {
|
| 2911 |
2960 |
blossoms.push_back(c);
|
| 2912 |
2961 |
}
|
| 2913 |
2962 |
|
| 2914 |
2963 |
for (int i = 0; i < int(blossoms.size()); ++i) {
|
| 2915 |
2964 |
|
| 2916 |
2965 |
Value offset = (*_blossom_data)[blossoms[i]].offset;
|
| 2917 |
2966 |
(*_blossom_data)[blossoms[i]].pot += 2 * offset;
|
| 2918 |
2967 |
for (typename BlossomSet::ItemIt n(*_blossom_set, blossoms[i]);
|
| 2919 |
2968 |
n != INVALID; ++n) {
|
| 2920 |
2969 |
(*_node_data)[(*_node_index)[n]].pot -= offset;
|
| 2921 |
2970 |
}
|
| 2922 |
2971 |
|
| 2923 |
2972 |
Arc matching = (*_blossom_data)[blossoms[i]].next;
|
| 2924 |
2973 |
Node base = _graph.source(matching);
|
| 2925 |
2974 |
extractBlossom(blossoms[i], base, matching);
|
| 2926 |
2975 |
}
|
| 2927 |
2976 |
}
|
| 2928 |
2977 |
|
| 2929 |
2978 |
public:
|
| 2930 |
2979 |
|
| 2931 |
2980 |
/// \brief Constructor
|
| 2932 |
2981 |
///
|
| 2933 |
2982 |
/// Constructor.
|
| 2934 |
2983 |
MaxWeightedPerfectMatching(const Graph& graph, const WeightMap& weight)
|
| 2935 |
2984 |
: _graph(graph), _weight(weight), _matching(0),
|
| 2936 |
2985 |
_node_potential(0), _blossom_potential(), _blossom_node_list(),
|
| 2937 |
2986 |
_node_num(0), _blossom_num(0),
|
| 2938 |
2987 |
|
| 2939 |
2988 |
_blossom_index(0), _blossom_set(0), _blossom_data(0),
|
| 2940 |
2989 |
_node_index(0), _node_heap_index(0), _node_data(0),
|
| 2941 |
2990 |
_tree_set_index(0), _tree_set(0),
|
| 2942 |
2991 |
|
| 2943 |
2992 |
_delta2_index(0), _delta2(0),
|
| 2944 |
2993 |
_delta3_index(0), _delta3(0),
|
| 2945 |
2994 |
_delta4_index(0), _delta4(0),
|
| 2946 |
2995 |
|
| 2947 |
2996 |
_delta_sum(), _unmatched(0),
|
| 2948 |
2997 |
|
| 2949 |
2998 |
_fractional(0)
|
| 2950 |
2999 |
{}
|
| 2951 |
3000 |
|
| 2952 |
3001 |
~MaxWeightedPerfectMatching() {
|
| 2953 |
3002 |
destroyStructures();
|
| 2954 |
3003 |
if (_fractional) {
|
| 2955 |
3004 |
delete _fractional;
|
| 2956 |
3005 |
}
|
| 2957 |
3006 |
}
|
| 2958 |
3007 |
|
| 2959 |
3008 |
/// \name Execution Control
|
| 2960 |
3009 |
/// The simplest way to execute the algorithm is to use the
|
| 2961 |
3010 |
/// \ref run() member function.
|
| 2962 |
3011 |
|
| 2963 |
3012 |
///@{
|
| 2964 |
3013 |
|
| 2965 |
3014 |
/// \brief Initialize the algorithm
|
| 2966 |
3015 |
///
|
| 2967 |
3016 |
/// This function initializes the algorithm.
|
| 2968 |
3017 |
void init() {
|
| 2969 |
3018 |
createStructures();
|
| 2970 |
3019 |
|
|
3020 |
_blossom_node_list.clear();
|
|
3021 |
_blossom_potential.clear();
|
|
3022 |
|
| 2971 |
3023 |
for (ArcIt e(_graph); e != INVALID; ++e) {
|
| 2972 |
3024 |
(*_node_heap_index)[e] = BinHeap<Value, IntArcMap>::PRE_HEAP;
|
| 2973 |
3025 |
}
|
| 2974 |
3026 |
for (EdgeIt e(_graph); e != INVALID; ++e) {
|
| 2975 |
3027 |
(*_delta3_index)[e] = _delta3->PRE_HEAP;
|
| 2976 |
3028 |
}
|
| 2977 |
3029 |
for (int i = 0; i < _blossom_num; ++i) {
|
| 2978 |
3030 |
(*_delta2_index)[i] = _delta2->PRE_HEAP;
|
| 2979 |
3031 |
(*_delta4_index)[i] = _delta4->PRE_HEAP;
|
| 2980 |
3032 |
}
|
| 2981 |
3033 |
|
| 2982 |
3034 |
_unmatched = _node_num;
|
| 2983 |
3035 |
|
|
3036 |
_delta2->clear();
|
|
3037 |
_delta3->clear();
|
|
3038 |
_delta4->clear();
|
|
3039 |
_blossom_set->clear();
|
|
3040 |
_tree_set->clear();
|
|
3041 |
|
| 2984 |
3042 |
int index = 0;
|
| 2985 |
3043 |
for (NodeIt n(_graph); n != INVALID; ++n) {
|
| 2986 |
3044 |
Value max = - std::numeric_limits<Value>::max();
|
| 2987 |
3045 |
for (OutArcIt e(_graph, n); e != INVALID; ++e) {
|
| 2988 |
3046 |
if (_graph.target(e) == n) continue;
|
| 2989 |
3047 |
if ((dualScale * _weight[e]) / 2 > max) {
|
| 2990 |
3048 |
max = (dualScale * _weight[e]) / 2;
|
| 2991 |
3049 |
}
|
| 2992 |
3050 |
}
|
| 2993 |
3051 |
(*_node_index)[n] = index;
|
|
3052 |
(*_node_data)[index].heap_index.clear();
|
|
3053 |
(*_node_data)[index].heap.clear();
|
| 2994 |
3054 |
(*_node_data)[index].pot = max;
|
| 2995 |
3055 |
int blossom =
|
| 2996 |
3056 |
_blossom_set->insert(n, std::numeric_limits<Value>::max());
|
| 2997 |
3057 |
|
| 2998 |
3058 |
_tree_set->insert(blossom);
|
| 2999 |
3059 |
|
| 3000 |
3060 |
(*_blossom_data)[blossom].status = EVEN;
|
| 3001 |
3061 |
(*_blossom_data)[blossom].pred = INVALID;
|
| 3002 |
3062 |
(*_blossom_data)[blossom].next = INVALID;
|
| 3003 |
3063 |
(*_blossom_data)[blossom].pot = 0;
|
| 3004 |
3064 |
(*_blossom_data)[blossom].offset = 0;
|
| 3005 |
3065 |
++index;
|
| 3006 |
3066 |
}
|
| 3007 |
3067 |
for (EdgeIt e(_graph); e != INVALID; ++e) {
|
| 3008 |
3068 |
int si = (*_node_index)[_graph.u(e)];
|
| 3009 |
3069 |
int ti = (*_node_index)[_graph.v(e)];
|
| 3010 |
3070 |
if (_graph.u(e) != _graph.v(e)) {
|
| 3011 |
3071 |
_delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot -
|
| 3012 |
3072 |
dualScale * _weight[e]) / 2);
|
| 3013 |
3073 |
}
|
| 3014 |
3074 |
}
|
| 3015 |
3075 |
}
|
| 3016 |
3076 |
|
| 3017 |
3077 |
/// \brief Initialize the algorithm with fractional matching
|
| 3018 |
3078 |
///
|
| 3019 |
3079 |
/// This function initializes the algorithm with a fractional
|
| 3020 |
3080 |
/// matching. This initialization is also called jumpstart heuristic.
|
| 3021 |
3081 |
void fractionalInit() {
|
| 3022 |
3082 |
createStructures();
|
| 3023 |
3083 |
|
| 3024 |
3084 |
if (_fractional == 0) {
|
| 3025 |
3085 |
_fractional = new FractionalMatching(_graph, _weight, false);
|
| 3026 |
3086 |
}
|
| 3027 |
3087 |
if (!_fractional->run()) {
|
| 3028 |
3088 |
_unmatched = -1;
|
| 3029 |
3089 |
return;
|
| 3030 |
3090 |
}
|
| 3031 |
3091 |
|
| 3032 |
3092 |
for (ArcIt e(_graph); e != INVALID; ++e) {
|
| 3033 |
3093 |
(*_node_heap_index)[e] = BinHeap<Value, IntArcMap>::PRE_HEAP;
|
| 3034 |
3094 |
}
|
| 3035 |
3095 |
for (EdgeIt e(_graph); e != INVALID; ++e) {
|
| 3036 |
3096 |
(*_delta3_index)[e] = _delta3->PRE_HEAP;
|
| 3037 |
3097 |
}
|
| 3038 |
3098 |
for (int i = 0; i < _blossom_num; ++i) {
|
| 3039 |
3099 |
(*_delta2_index)[i] = _delta2->PRE_HEAP;
|
| 3040 |
3100 |
(*_delta4_index)[i] = _delta4->PRE_HEAP;
|
| 3041 |
3101 |
}
|
| 3042 |
3102 |
|
| 3043 |
3103 |
_unmatched = 0;
|
| 3044 |
3104 |
|
| 3045 |
3105 |
int index = 0;
|
| 3046 |
3106 |
for (NodeIt n(_graph); n != INVALID; ++n) {
|
| 3047 |
3107 |
Value pot = _fractional->nodeValue(n);
|
| 3048 |
3108 |
(*_node_index)[n] = index;
|
| 3049 |
3109 |
(*_node_data)[index].pot = pot;
|
| 3050 |
3110 |
int blossom =
|
| 3051 |
3111 |
_blossom_set->insert(n, std::numeric_limits<Value>::max());
|
| 3052 |
3112 |
|
| 3053 |
3113 |
(*_blossom_data)[blossom].status = MATCHED;
|
| 3054 |
3114 |
(*_blossom_data)[blossom].pred = INVALID;
|
| 3055 |
3115 |
(*_blossom_data)[blossom].next = _fractional->matching(n);
|
| 3056 |
3116 |
(*_blossom_data)[blossom].pot = 0;
|
| 3057 |
3117 |
(*_blossom_data)[blossom].offset = 0;
|
| 3058 |
3118 |
++index;
|
| 3059 |
3119 |
}
|
| 3060 |
3120 |
|
| 3061 |
3121 |
typename Graph::template NodeMap<bool> processed(_graph, false);
|
| 3062 |
3122 |
for (NodeIt n(_graph); n != INVALID; ++n) {
|
| 3063 |
3123 |
if (processed[n]) continue;
|
| 3064 |
3124 |
processed[n] = true;
|
| 3065 |
3125 |
if (_fractional->matching(n) == INVALID) continue;
|
| 3066 |
3126 |
int num = 1;
|
| 3067 |
3127 |
Node v = _graph.target(_fractional->matching(n));
|
| 3068 |
3128 |
while (n != v) {
|
| 3069 |
3129 |
processed[v] = true;
|
| 3070 |
3130 |
v = _graph.target(_fractional->matching(v));
|
| 3071 |
3131 |
++num;
|
| 3072 |
3132 |
}
|
| 3073 |
3133 |
|
| 3074 |
3134 |
if (num % 2 == 1) {
|
| 3075 |
3135 |
std::vector<int> subblossoms(num);
|
| 3076 |
3136 |
|
| 3077 |
3137 |
subblossoms[--num] = _blossom_set->find(n);
|
| 3078 |
3138 |
v = _graph.target(_fractional->matching(n));
|
| 3079 |
3139 |
while (n != v) {
|
| 3080 |
3140 |
subblossoms[--num] = _blossom_set->find(v);
|
| 3081 |
3141 |
v = _graph.target(_fractional->matching(v));
|
| 3082 |
3142 |
}
|
| 3083 |
3143 |
|
| 3084 |
3144 |
int surface =
|
| 3085 |
3145 |
_blossom_set->join(subblossoms.begin(), subblossoms.end());
|
| 3086 |
3146 |
(*_blossom_data)[surface].status = EVEN;
|
| 3087 |
3147 |
(*_blossom_data)[surface].pred = INVALID;
|
| 3088 |
3148 |
(*_blossom_data)[surface].next = INVALID;
|
| 3089 |
3149 |
(*_blossom_data)[surface].pot = 0;
|
| 3090 |
3150 |
(*_blossom_data)[surface].offset = 0;
|
| 3091 |
3151 |
|
| 3092 |
3152 |
_tree_set->insert(surface);
|
| 3093 |
3153 |
++_unmatched;
|
| 3094 |
3154 |
}
|
| 3095 |
3155 |
}
|
| 3096 |
3156 |
|
| 3097 |
3157 |
for (EdgeIt e(_graph); e != INVALID; ++e) {
|
| 3098 |
3158 |
int si = (*_node_index)[_graph.u(e)];
|
| 3099 |
3159 |
int sb = _blossom_set->find(_graph.u(e));
|
| 3100 |
3160 |
int ti = (*_node_index)[_graph.v(e)];
|
| 3101 |
3161 |
int tb = _blossom_set->find(_graph.v(e));
|
| 3102 |
3162 |
if ((*_blossom_data)[sb].status == EVEN &&
|
| 3103 |
3163 |
(*_blossom_data)[tb].status == EVEN && sb != tb) {
|
| 3104 |
3164 |
_delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot -
|
| 3105 |
3165 |
dualScale * _weight[e]) / 2);
|
| 3106 |
3166 |
}
|
| 3107 |
3167 |
}
|
| 3108 |
3168 |
|
| 3109 |
3169 |
for (NodeIt n(_graph); n != INVALID; ++n) {
|
| 3110 |
3170 |
int nb = _blossom_set->find(n);
|
| 3111 |
3171 |
if ((*_blossom_data)[nb].status != MATCHED) continue;
|
| 3112 |
3172 |
int ni = (*_node_index)[n];
|
| 3113 |
3173 |
|
| 3114 |
3174 |
for (OutArcIt e(_graph, n); e != INVALID; ++e) {
|
| 3115 |
3175 |
Node v = _graph.target(e);
|
| 3116 |
3176 |
int vb = _blossom_set->find(v);
|
| 3117 |
3177 |
int vi = (*_node_index)[v];
|
| 3118 |
3178 |
|
| 3119 |
3179 |
Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
|
| 3120 |
3180 |
dualScale * _weight[e];
|
| 3121 |
3181 |
|
| 3122 |
3182 |
if ((*_blossom_data)[vb].status == EVEN) {
|
| 3123 |
3183 |
|
| 3124 |
3184 |
int vt = _tree_set->find(vb);
|
| 3125 |
3185 |
|
| 3126 |
3186 |
typename std::map<int, Arc>::iterator it =
|
| 3127 |
3187 |
(*_node_data)[ni].heap_index.find(vt);
|
| 3128 |
3188 |
|
| 3129 |
3189 |
if (it != (*_node_data)[ni].heap_index.end()) {
|
| 3130 |
3190 |
if ((*_node_data)[ni].heap[it->second] > rw) {
|
| 3131 |
3191 |
(*_node_data)[ni].heap.replace(it->second, e);
|
| 3132 |
3192 |
(*_node_data)[ni].heap.decrease(e, rw);
|
| 3133 |
3193 |
it->second = e;
|
| 3134 |
3194 |
}
|
| 3135 |
3195 |
} else {
|
| 3136 |
3196 |
(*_node_data)[ni].heap.push(e, rw);
|
| 3137 |
3197 |
(*_node_data)[ni].heap_index.insert(std::make_pair(vt, e));
|
| 3138 |
3198 |
}
|
| 3139 |
3199 |
}
|
| 3140 |
3200 |
}
|
| 3141 |
3201 |
|
| 3142 |
3202 |
if (!(*_node_data)[ni].heap.empty()) {
|
| 3143 |
3203 |
_blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
|
| 3144 |
3204 |
_delta2->push(nb, _blossom_set->classPrio(nb));
|
| 3145 |
3205 |
}
|
| 3146 |
3206 |
}
|
| 3147 |
3207 |
}
|
| 3148 |
3208 |
|
| 3149 |
3209 |
/// \brief Start the algorithm
|
| 3150 |
3210 |
///
|
| 3151 |
3211 |
/// This function starts the algorithm.
|
| 3152 |
3212 |
///
|
| 3153 |
3213 |
/// \pre \ref init() or \ref fractionalInit() must be called before
|
| 3154 |
3214 |
/// using this function.
|
| 3155 |
3215 |
bool start() {
|
| 3156 |
3216 |
enum OpType {
|
| 3157 |
3217 |
D2, D3, D4
|
| 3158 |
3218 |
};
|
| 3159 |
3219 |
|
| 3160 |
3220 |
if (_unmatched == -1) return false;
|
| 3161 |
3221 |
|
| 3162 |
3222 |
while (_unmatched > 0) {
|
| 3163 |
3223 |
Value d2 = !_delta2->empty() ?
|
| 3164 |
3224 |
_delta2->prio() : std::numeric_limits<Value>::max();
|
| 3165 |
3225 |
|
| 3166 |
3226 |
Value d3 = !_delta3->empty() ?
|
| 3167 |
3227 |
_delta3->prio() : std::numeric_limits<Value>::max();
|
| 3168 |
3228 |
|
| 3169 |
3229 |
Value d4 = !_delta4->empty() ?
|
| 3170 |
3230 |
_delta4->prio() : std::numeric_limits<Value>::max();
|
| 3171 |
3231 |
|
| 3172 |
3232 |
_delta_sum = d3; OpType ot = D3;
|
| 3173 |
3233 |
if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; }
|
| 3174 |
3234 |
if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; }
|
| 3175 |
3235 |
|
| 3176 |
3236 |
if (_delta_sum == std::numeric_limits<Value>::max()) {
|
| 3177 |
3237 |
return false;
|
| 3178 |
3238 |
}
|
| 3179 |
3239 |
|
| 3180 |
3240 |
switch (ot) {
|
| 3181 |
3241 |
case D2:
|
| 3182 |
3242 |
{
|
| 3183 |
3243 |
int blossom = _delta2->top();
|
| 3184 |
3244 |
Node n = _blossom_set->classTop(blossom);
|
| 3185 |
3245 |
Arc e = (*_node_data)[(*_node_index)[n]].heap.top();
|
| 3186 |
3246 |
extendOnArc(e);
|
| 3187 |
3247 |
}
|
| 3188 |
3248 |
break;
|
| 3189 |
3249 |
case D3:
|
| 3190 |
3250 |
{
|
| 3191 |
3251 |
Edge e = _delta3->top();
|
| 3192 |
3252 |
|
| 3193 |
3253 |
int left_blossom = _blossom_set->find(_graph.u(e));
|
| 3194 |
3254 |
int right_blossom = _blossom_set->find(_graph.v(e));
|
| 3195 |
3255 |
|
| 3196 |
3256 |
if (left_blossom == right_blossom) {
|
| 3197 |
3257 |
_delta3->pop();
|
| 3198 |
3258 |
} else {
|
| 3199 |
3259 |
int left_tree = _tree_set->find(left_blossom);
|
| 3200 |
3260 |
int right_tree = _tree_set->find(right_blossom);
|
| 3201 |
3261 |
|
| 3202 |
3262 |
if (left_tree == right_tree) {
|
| 3203 |
3263 |
shrinkOnEdge(e, left_tree);
|
| 3204 |
3264 |
} else {
|
| 3205 |
3265 |
augmentOnEdge(e);
|
| 3206 |
3266 |
_unmatched -= 2;
|
| 3207 |
3267 |
}
|
| 3208 |
3268 |
}
|
| 3209 |
3269 |
} break;
|
| 3210 |
3270 |
case D4:
|
| 3211 |
3271 |
splitBlossom(_delta4->top());
|
| 3212 |
3272 |
break;
|
| 3213 |
3273 |
}
|
| 3214 |
3274 |
}
|
| 3215 |
3275 |
extractMatching();
|
| 3216 |
3276 |
return true;
|
| 3217 |
3277 |
}
|
| 3218 |
3278 |
|
| 3219 |
3279 |
/// \brief Run the algorithm.
|
| 3220 |
3280 |
///
|
| 3221 |
3281 |
/// This method runs the \c %MaxWeightedPerfectMatching algorithm.
|
| 3222 |
3282 |
///
|
| 3223 |
3283 |
/// \note mwpm.run() is just a shortcut of the following code.
|
| 3224 |
3284 |
/// \code
|
| 3225 |
3285 |
/// mwpm.fractionalInit();
|
| 3226 |
3286 |
/// mwpm.start();
|
| 3227 |
3287 |
/// \endcode
|
| 3228 |
3288 |
bool run() {
|
| 3229 |
3289 |
fractionalInit();
|
| 3230 |
3290 |
return start();
|
| 3231 |
3291 |
}
|
| 3232 |
3292 |
|
| 3233 |
3293 |
/// @}
|
| 3234 |
3294 |
|
| 3235 |
3295 |
/// \name Primal Solution
|
| 3236 |
3296 |
/// Functions to get the primal solution, i.e. the maximum weighted
|
| 3237 |
3297 |
/// perfect matching.\n
|
| 3238 |
3298 |
/// Either \ref run() or \ref start() function should be called before
|
| 3239 |
3299 |
/// using them.
|
| 3240 |
3300 |
|
| 3241 |
3301 |
/// @{
|
| 3242 |
3302 |
|
| 3243 |
3303 |
/// \brief Return the weight of the matching.
|
| 3244 |
3304 |
///
|
| 3245 |
3305 |
/// This function returns the weight of the found matching.
|
| 3246 |
3306 |
///
|
| 3247 |
3307 |
/// \pre Either run() or start() must be called before using this function.
|
| 3248 |
3308 |
Value matchingWeight() const {
|
| 3249 |
3309 |
Value sum = 0;
|
| 3250 |
3310 |
for (NodeIt n(_graph); n != INVALID; ++n) {
|
| 3251 |
3311 |
if ((*_matching)[n] != INVALID) {
|
| 3252 |
3312 |
sum += _weight[(*_matching)[n]];
|
| 3253 |
3313 |
}
|
| 3254 |
3314 |
}
|
| 3255 |
3315 |
return sum / 2;
|
| 3256 |
3316 |
}
|
| 3257 |
3317 |
|
| 3258 |
3318 |
/// \brief Return \c true if the given edge is in the matching.
|
| 3259 |
3319 |
///
|
| 3260 |
3320 |
/// This function returns \c true if the given edge is in the found
|
| 3261 |
3321 |
/// matching.
|
| 3262 |
3322 |
///
|
| 3263 |
3323 |
/// \pre Either run() or start() must be called before using this function.
|
| 3264 |
3324 |
bool matching(const Edge& edge) const {
|
| 3265 |
3325 |
return static_cast<const Edge&>((*_matching)[_graph.u(edge)]) == edge;
|
| 3266 |
3326 |
}
|
| 3267 |
3327 |
|
| 3268 |
3328 |
/// \brief Return the matching arc (or edge) incident to the given node.
|
| 3269 |
3329 |
///
|
| 3270 |
3330 |
/// This function returns the matching arc (or edge) incident to the
|
| 3271 |
3331 |
/// given node in the found matching or \c INVALID if the node is
|
| 3272 |
3332 |
/// not covered by the matching.
|
| 3273 |
3333 |
///
|
| 3274 |
3334 |
/// \pre Either run() or start() must be called before using this function.
|
| 3275 |
3335 |
Arc matching(const Node& node) const {
|
| 3276 |
3336 |
return (*_matching)[node];
|
| 3277 |
3337 |
}
|
| 3278 |
3338 |
|
| 3279 |
3339 |
/// \brief Return a const reference to the matching map.
|
| 3280 |
3340 |
///
|
| 3281 |
3341 |
/// This function returns a const reference to a node map that stores
|
| 3282 |
3342 |
/// the matching arc (or edge) incident to each node.
|
| 3283 |
3343 |
const MatchingMap& matchingMap() const {
|
| 3284 |
3344 |
return *_matching;
|
| 3285 |
3345 |
}
|
| 3286 |
3346 |
|
| 3287 |
3347 |
/// \brief Return the mate of the given node.
|
| 3288 |
3348 |
///
|
| 3289 |
3349 |
/// This function returns the mate of the given node in the found
|
| 3290 |
3350 |
/// matching or \c INVALID if the node is not covered by the matching.
|
| 3291 |
3351 |
///
|
| 3292 |
3352 |
/// \pre Either run() or start() must be called before using this function.
|
| 3293 |
3353 |
Node mate(const Node& node) const {
|
| 3294 |
3354 |
return _graph.target((*_matching)[node]);
|
| 3295 |
3355 |
}
|
| 3296 |
3356 |
|
| 3297 |
3357 |
/// @}
|
| 3298 |
3358 |
|
| 3299 |
3359 |
/// \name Dual Solution
|
| 3300 |
3360 |
/// Functions to get the dual solution.\n
|
| 3301 |
3361 |
/// Either \ref run() or \ref start() function should be called before
|
| 3302 |
3362 |
/// using them.
|
| 3303 |
3363 |
|
| 3304 |
3364 |
/// @{
|
| 3305 |
3365 |
|
| 3306 |
3366 |
/// \brief Return the value of the dual solution.
|
| 3307 |
3367 |
///
|
| 3308 |
3368 |
/// This function returns the value of the dual solution.
|
| 3309 |
3369 |
/// It should be equal to the primal value scaled by \ref dualScale
|
| 3310 |
3370 |
/// "dual scale".
|
| 3311 |
3371 |
///
|
| 3312 |
3372 |
/// \pre Either run() or start() must be called before using this function.
|
| 3313 |
3373 |
Value dualValue() const {
|
| 3314 |
3374 |
Value sum = 0;
|
| 3315 |
3375 |
for (NodeIt n(_graph); n != INVALID; ++n) {
|
| 3316 |
3376 |
sum += nodeValue(n);
|
| 3317 |
3377 |
}
|
| 3318 |
3378 |
for (int i = 0; i < blossomNum(); ++i) {
|
| 3319 |
3379 |
sum += blossomValue(i) * (blossomSize(i) / 2);
|
| 3320 |
3380 |
}
|
| 3321 |
3381 |
return sum;
|
| 3322 |
3382 |
}
|
| 3323 |
3383 |
|
| 3324 |
3384 |
/// \brief Return the dual value (potential) of the given node.
|
| 3325 |
3385 |
///
|
| 3326 |
3386 |
/// This function returns the dual value (potential) of the given node.
|
| 3327 |
3387 |
///
|
| 3328 |
3388 |
/// \pre Either run() or start() must be called before using this function.
|
| 3329 |
3389 |
Value nodeValue(const Node& n) const {
|
| 3330 |
3390 |
return (*_node_potential)[n];
|
| 3331 |
3391 |
}
|
| 3332 |
3392 |
|
| 3333 |
3393 |
/// \brief Return the number of the blossoms in the basis.
|
| 3334 |
3394 |
///
|
| 3335 |
3395 |
/// This function returns the number of the blossoms in the basis.
|
| 3336 |
3396 |
///
|
| 3337 |
3397 |
/// \pre Either run() or start() must be called before using this function.
|
| 3338 |
3398 |
/// \see BlossomIt
|
| 3339 |
3399 |
int blossomNum() const {
|
| 3340 |
3400 |
return _blossom_potential.size();
|
| 3341 |
3401 |
}
|
| 3342 |
3402 |
|
| 3343 |
3403 |
/// \brief Return the number of the nodes in the given blossom.
|
| 3344 |
3404 |
///
|
| 3345 |
3405 |
/// This function returns the number of the nodes in the given blossom.
|
| 3346 |
3406 |
///
|
| 3347 |
3407 |
/// \pre Either run() or start() must be called before using this function.
|
| 3348 |
3408 |
/// \see BlossomIt
|
| 3349 |
3409 |
int blossomSize(int k) const {
|
| 3350 |
3410 |
return _blossom_potential[k].end - _blossom_potential[k].begin;
|
| 3351 |
3411 |
}
|
| 3352 |
3412 |
|
| 3353 |
3413 |
/// \brief Return the dual value (ptential) of the given blossom.
|
| 3354 |
3414 |
///
|
| 3355 |
3415 |
/// This function returns the dual value (ptential) of the given blossom.
|
| 3356 |
3416 |
///
|
| 3357 |
3417 |
/// \pre Either run() or start() must be called before using this function.
|
| 3358 |
3418 |
Value blossomValue(int k) const {
|
| 3359 |
3419 |
return _blossom_potential[k].value;
|
| 3360 |
3420 |
}
|
| 3361 |
3421 |
|
| 3362 |
3422 |
/// \brief Iterator for obtaining the nodes of a blossom.
|
| 3363 |
3423 |
///
|
| 3364 |
3424 |
/// This class provides an iterator for obtaining the nodes of the
|
| 3365 |
3425 |
/// given blossom. It lists a subset of the nodes.
|
| 3366 |
3426 |
/// Before using this iterator, you must allocate a
|
| 3367 |
3427 |
/// MaxWeightedPerfectMatching class and execute it.
|
| 3368 |
3428 |
class BlossomIt {
|
| 3369 |
3429 |
public:
|
| 3370 |
3430 |
|
| 3371 |
3431 |
/// \brief Constructor.
|
| 3372 |
3432 |
///
|
| 3373 |
3433 |
/// Constructor to get the nodes of the given variable.
|
| 3374 |
3434 |
///
|
| 3375 |
3435 |
/// \pre Either \ref MaxWeightedPerfectMatching::run() "algorithm.run()"
|
| 3376 |
3436 |
/// or \ref MaxWeightedPerfectMatching::start() "algorithm.start()"
|
| 3377 |
3437 |
/// must be called before initializing this iterator.
|