Changeset 2051:08652c1763f6 in lemon-0.x
- Timestamp:
- 04/14/06 20:07:33 (19 years ago)
- Branch:
- default
- Phase:
- public
- Convert:
- svn:c9d7d8f5-90d6-0310-b91f-818b3a526b0e/lemon/trunk@2694
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
lemon/bipartite_matching.h
r2040 r2051 20 20 #define LEMON_BIPARTITE_MATCHING 21 21 22 #include <lemon/bpugraph_adaptor.h> 23 #include <lemon/bfs.h> 22 #include <functional> 23 24 #include <lemon/bin_heap.h> 25 #include <lemon/maps.h> 24 26 25 27 #include <iostream> … … 36 38 /// 37 39 /// Bipartite Max Cardinality Matching algorithm. This class implements 38 /// the Hopcroft-Karp algorithm w ich has \f$ O(e\sqrt{n}) \f$ time40 /// the Hopcroft-Karp algorithm which has \f$ O(e\sqrt{n}) \f$ time 39 41 /// complexity. 40 42 template <typename BpUGraph> … … 84 86 bnode_matching[it] = INVALID; 85 87 } 86 matching_ value = 0;88 matching_size = 0; 87 89 } 88 90 … … 93 95 /// the matching than from the initial empty matching. 94 96 void greedyInit() { 95 matching_ value = 0;97 matching_size = 0; 96 98 for (BNodeIt it(*graph); it != INVALID; ++it) { 97 99 bnode_matching[it] = INVALID; … … 103 105 anode_matching[it] = jt; 104 106 bnode_matching[graph->bNode(jt)] = jt; 105 ++matching_ value;107 ++matching_size; 106 108 break; 107 109 } … … 121 123 bnode_matching[it] = INVALID; 122 124 } 123 matching_ value = 0;125 matching_size = 0; 124 126 for (UEdgeIt it(*graph); it != INVALID; ++it) { 125 127 if (matching[it]) { 126 ++matching_ value;128 ++matching_size; 127 129 anode_matching[graph->aNode(it)] = it; 128 130 bnode_matching[graph->bNode(it)] = it; … … 143 145 bnode_matching[it] = INVALID; 144 146 } 145 matching_ value = 0;147 matching_size = 0; 146 148 for (UEdgeIt it(*graph); it != INVALID; ++it) { 147 149 if (matching[it]) { 148 ++matching_ value;150 ++matching_size; 149 151 if (anode_matching[graph->aNode(it)] != INVALID) { 150 152 return false; … … 247 249 aused[anode] = true; 248 250 } 249 ++matching_ value;251 ++matching_size; 250 252 251 253 } … … 298 300 299 301 } 300 ++matching_ value;302 ++matching_size; 301 303 return true; 302 304 } else { … … 408 410 } 409 411 } 410 return matching_ value;412 return matching_size; 411 413 } 412 414 … … 420 422 matching[it] = it == anode_matching[graph->aNode(it)]; 421 423 } 422 return matching_ value;424 return matching_size; 423 425 } 424 426 … … 446 448 /// 447 449 /// Gives back the number of the matching edges. 448 int matching Value() const {449 return matching_ value;450 int matchingSize() const { 451 return matching_size; 450 452 } 451 453 … … 458 460 const Graph *graph; 459 461 460 int matching_ value;462 int matching_size; 461 463 462 464 }; 463 465 466 /// \brief Default traits class for weighted bipartite matching algoritms. 467 /// 468 /// Default traits class for weighted bipartite matching algoritms. 469 /// \param _BpUGraph The bipartite undirected graph type. 470 /// \param _WeightMap Type of weight map. 471 template <typename _BpUGraph, typename _WeightMap> 472 struct WeightedBipartiteMatchingDefaultTraits { 473 /// \brief The type of the weight of the undirected edges. 474 typedef typename _WeightMap::Value Value; 475 476 /// The undirected bipartite graph type the algorithm runs on. 477 typedef _BpUGraph BpUGraph; 478 479 /// The map of the edges weights 480 typedef _WeightMap WeightMap; 481 482 /// \brief The cross reference type used by heap. 483 /// 484 /// The cross reference type used by heap. 485 /// Usually it is \c Graph::NodeMap<int>. 486 typedef typename BpUGraph::template NodeMap<int> HeapCrossRef; 487 488 /// \brief Instantiates a HeapCrossRef. 489 /// 490 /// This function instantiates a \ref HeapCrossRef. 491 /// \param graph is the graph, to which we would like to define the 492 /// HeapCrossRef. 493 static HeapCrossRef *createHeapCrossRef(const BpUGraph &graph) { 494 return new HeapCrossRef(graph); 495 } 496 497 /// \brief The heap type used by weighted matching algorithms. 498 /// 499 /// The heap type used by weighted matching algorithms. It should 500 /// minimize the priorities and the heap's key type is the graph's 501 /// anode graph's node. 502 /// 503 /// \sa BinHeap 504 typedef BinHeap<typename BpUGraph::Node, Value, HeapCrossRef> Heap; 505 506 /// \brief Instantiates a Heap. 507 /// 508 /// This function instantiates a \ref Heap. 509 /// \param crossref The cross reference of the heap. 510 static Heap *createHeap(HeapCrossRef& crossref) { 511 return new Heap(crossref); 512 } 513 514 }; 515 516 517 /// \ingroup matching 518 /// 519 /// \brief Bipartite Max Weighted Matching algorithm 520 /// 521 /// This class implements the bipartite Max Weighted Matching 522 /// algorithm. It uses the successive shortest path algorithm to 523 /// calculate the maximum weighted matching in the bipartite 524 /// graph. The algorithm can be used also to calculate the maximum 525 /// cardinality maximum weighted matching. The time complexity 526 /// of the algorithm is \f$ O(ne\log(n)) \f$ with the default binary 527 /// heap implementation but this can be improved to 528 /// \f$ O(n^2\log(n)+ne) \f$ if we use fibonacci heaps. 529 /// 530 /// The algorithm also provides a potential function on the nodes 531 /// which a dual solution of the matching algorithm and it can be 532 /// used to proof the optimality of the given pimal solution. 533 #ifdef DOXYGEN 534 template <typename _BpUGraph, typename _WeightMap, typename _Traits> 535 #else 536 template <typename _BpUGraph, 537 typename _WeightMap = typename _BpUGraph::template UEdgeMap<int>, 538 typename _Traits = WeightedBipartiteMatchingDefaultTraits<_BpUGraph, _WeightMap> > 539 #endif 540 class MaxWeightedBipartiteMatching { 541 public: 542 543 typedef _Traits Traits; 544 typedef typename Traits::BpUGraph BpUGraph; 545 typedef typename Traits::WeightMap WeightMap; 546 typedef typename Traits::Value Value; 547 548 protected: 549 550 typedef typename Traits::HeapCrossRef HeapCrossRef; 551 typedef typename Traits::Heap Heap; 552 553 554 typedef typename BpUGraph::Node Node; 555 typedef typename BpUGraph::ANodeIt ANodeIt; 556 typedef typename BpUGraph::BNodeIt BNodeIt; 557 typedef typename BpUGraph::UEdge UEdge; 558 typedef typename BpUGraph::UEdgeIt UEdgeIt; 559 typedef typename BpUGraph::IncEdgeIt IncEdgeIt; 560 561 typedef typename BpUGraph::template ANodeMap<UEdge> ANodeMatchingMap; 562 typedef typename BpUGraph::template BNodeMap<UEdge> BNodeMatchingMap; 563 564 typedef typename BpUGraph::template ANodeMap<Value> ANodePotentialMap; 565 typedef typename BpUGraph::template BNodeMap<Value> BNodePotentialMap; 566 567 568 public: 569 570 /// \brief \ref Exception for uninitialized parameters. 571 /// 572 /// This error represents problems in the initialization 573 /// of the parameters of the algorithms. 574 class UninitializedParameter : public lemon::UninitializedParameter { 575 public: 576 virtual const char* exceptionName() const { 577 return "lemon::MaxWeightedBipartiteMatching::UninitializedParameter"; 578 } 579 }; 580 581 ///\name Named template parameters 582 583 ///@{ 584 585 template <class H, class CR> 586 struct DefHeapTraits : public Traits { 587 typedef CR HeapCrossRef; 588 typedef H Heap; 589 static HeapCrossRef *createHeapCrossRef(const BpUGraph &) { 590 throw UninitializedParameter(); 591 } 592 static Heap *createHeap(HeapCrossRef &) { 593 throw UninitializedParameter(); 594 } 595 }; 596 597 /// \brief \ref named-templ-param "Named parameter" for setting heap 598 /// and cross reference type 599 /// 600 /// \ref named-templ-param "Named parameter" for setting heap and cross 601 /// reference type 602 template <class H, class CR = typename BpUGraph::template NodeMap<int> > 603 struct DefHeap 604 : public MaxWeightedBipartiteMatching<BpUGraph, WeightMap, 605 DefHeapTraits<H, CR> > { 606 typedef MaxWeightedBipartiteMatching<BpUGraph, WeightMap, 607 DefHeapTraits<H, CR> > Create; 608 }; 609 610 template <class H, class CR> 611 struct DefStandardHeapTraits : public Traits { 612 typedef CR HeapCrossRef; 613 typedef H Heap; 614 static HeapCrossRef *createHeapCrossRef(const BpUGraph &graph) { 615 return new HeapCrossRef(graph); 616 } 617 static Heap *createHeap(HeapCrossRef &crossref) { 618 return new Heap(crossref); 619 } 620 }; 621 622 /// \brief \ref named-templ-param "Named parameter" for setting heap and 623 /// cross reference type with automatic allocation 624 /// 625 /// \ref named-templ-param "Named parameter" for setting heap and cross 626 /// reference type. It can allocate the heap and the cross reference 627 /// object if the cross reference's constructor waits for the graph as 628 /// parameter and the heap's constructor waits for the cross reference. 629 template <class H, class CR = typename BpUGraph::template NodeMap<int> > 630 struct DefStandardHeap 631 : public MaxWeightedBipartiteMatching<BpUGraph, WeightMap, 632 DefStandardHeapTraits<H, CR> > { 633 typedef MaxWeightedBipartiteMatching<BpUGraph, WeightMap, 634 DefStandardHeapTraits<H, CR> > 635 Create; 636 }; 637 638 ///@} 639 640 641 /// \brief Constructor. 642 /// 643 /// Constructor of the algorithm. 644 MaxWeightedBipartiteMatching(const BpUGraph& _graph, 645 const WeightMap& _weight) 646 : graph(&_graph), weight(&_weight), 647 anode_matching(_graph), bnode_matching(_graph), 648 anode_potential(_graph), bnode_potential(_graph), 649 _heap_cross_ref(0), local_heap_cross_ref(false), 650 _heap(0), local_heap(0) {} 651 652 /// \brief Destructor. 653 /// 654 /// Destructor of the algorithm. 655 ~MaxWeightedBipartiteMatching() { 656 destroyStructures(); 657 } 658 659 /// \brief Sets the heap and the cross reference used by algorithm. 660 /// 661 /// Sets the heap and the cross reference used by algorithm. 662 /// If you don't use this function before calling \ref run(), 663 /// it will allocate one. The destuctor deallocates this 664 /// automatically allocated map, of course. 665 /// \return \c (*this) 666 MaxWeightedBipartiteMatching& heap(Heap& heap, HeapCrossRef &crossRef) { 667 if(local_heap_cross_ref) { 668 delete _heap_cross_ref; 669 local_heap_cross_ref = false; 670 } 671 _heap_cross_ref = &crossRef; 672 if(local_heap) { 673 delete _heap; 674 local_heap = false; 675 } 676 _heap = &heap; 677 return *this; 678 } 679 680 /// \name Execution control 681 /// The simplest way to execute the algorithm is to use 682 /// one of the member functions called \c run(). 683 /// \n 684 /// If you need more control on the execution, 685 /// first you must call \ref init() or one alternative for it. 686 /// Finally \ref start() will perform the matching computation or 687 /// with step-by-step execution you can augment the solution. 688 689 /// @{ 690 691 /// \brief Initalize the data structures. 692 /// 693 /// It initalizes the data structures and creates an empty matching. 694 void init() { 695 initStructures(); 696 for (ANodeIt it(*graph); it != INVALID; ++it) { 697 anode_matching[it] = INVALID; 698 anode_potential[it] = 0; 699 } 700 for (BNodeIt it(*graph); it != INVALID; ++it) { 701 bnode_matching[it] = INVALID; 702 bnode_potential[it] = 0; 703 for (IncEdgeIt jt(*graph, it); jt != INVALID; ++jt) { 704 if (-(*weight)[jt] <= bnode_potential[it]) { 705 bnode_potential[it] = -(*weight)[jt]; 706 } 707 } 708 } 709 matching_value = 0; 710 matching_size = 0; 711 } 712 713 714 /// \brief An augmenting phase of the weighted matching algorithm 715 /// 716 /// It runs an augmenting phase of the weighted matching 717 /// algorithm. The phase finds the best augmenting path and 718 /// augments only on this paths. 719 /// 720 /// The algorithm consists at most 721 /// of \f$ O(n) \f$ phase and one phase is \f$ O(n\log(n)+e) \f$ 722 /// long with Fibonacci heap or \f$ O((n+e)\log(n)) \f$ long 723 /// with binary heap. 724 /// \param decrease If the given parameter true the matching value 725 /// can be decreased in the augmenting phase. If we would like 726 /// to calculate the maximum cardinality maximum weighted matching 727 /// then we should let the algorithm to decrease the matching 728 /// value in order to increase the number of the matching edges. 729 bool augment(bool decrease = false) { 730 731 typename BpUGraph::template BNodeMap<Value> bdist(*graph); 732 typename BpUGraph::template BNodeMap<UEdge> bpred(*graph, INVALID); 733 734 Node bestNode = INVALID; 735 Value bestValue = 0; 736 737 _heap->clear(); 738 for (ANodeIt it(*graph); it != INVALID; ++it) { 739 (*_heap_cross_ref)[it] = Heap::PRE_HEAP; 740 } 741 742 for (ANodeIt it(*graph); it != INVALID; ++it) { 743 if (anode_matching[it] == INVALID) { 744 _heap->push(it, 0); 745 } 746 } 747 748 Value bdistMax = 0; 749 while (!_heap->empty()) { 750 Node anode = _heap->top(); 751 Value avalue = _heap->prio(); 752 _heap->pop(); 753 for (IncEdgeIt jt(*graph, anode); jt != INVALID; ++jt) { 754 if (jt == anode_matching[anode]) continue; 755 Node bnode = graph->bNode(jt); 756 Value bvalue = avalue + anode_potential[anode] - 757 bnode_potential[bnode] - (*weight)[jt]; 758 if (bpred[bnode] == INVALID || bvalue < bdist[bnode]) { 759 bdist[bnode] = bvalue; 760 bpred[bnode] = jt; 761 } 762 if (bvalue > bdistMax) { 763 bdistMax = bvalue; 764 } 765 if (bnode_matching[bnode] != INVALID) { 766 Node newanode = graph->aNode(bnode_matching[bnode]); 767 switch (_heap->state(newanode)) { 768 case Heap::PRE_HEAP: 769 _heap->push(newanode, bvalue); 770 break; 771 case Heap::IN_HEAP: 772 if (bvalue < (*_heap)[newanode]) { 773 _heap->decrease(newanode, bvalue); 774 } 775 break; 776 case Heap::POST_HEAP: 777 break; 778 } 779 } else { 780 if (bestNode == INVALID || 781 - bvalue - bnode_potential[bnode] > bestValue) { 782 bestValue = - bvalue - bnode_potential[bnode]; 783 bestNode = bnode; 784 } 785 } 786 } 787 } 788 789 if (bestNode == INVALID || (!decrease && bestValue < 0)) { 790 return false; 791 } 792 793 matching_value += bestValue; 794 ++matching_size; 795 796 for (BNodeIt it(*graph); it != INVALID; ++it) { 797 if (bpred[it] != INVALID) { 798 bnode_potential[it] += bdist[it]; 799 } else { 800 bnode_potential[it] += bdistMax; 801 } 802 } 803 for (ANodeIt it(*graph); it != INVALID; ++it) { 804 if (anode_matching[it] != INVALID) { 805 Node bnode = graph->bNode(anode_matching[it]); 806 if (bpred[bnode] != INVALID) { 807 anode_potential[it] += bdist[bnode]; 808 } else { 809 anode_potential[it] += bdistMax; 810 } 811 } 812 } 813 814 while (bestNode != INVALID) { 815 UEdge uedge = bpred[bestNode]; 816 Node anode = graph->aNode(uedge); 817 818 bnode_matching[bestNode] = uedge; 819 if (anode_matching[anode] != INVALID) { 820 bestNode = graph->bNode(anode_matching[anode]); 821 } else { 822 bestNode = INVALID; 823 } 824 anode_matching[anode] = uedge; 825 } 826 827 828 return true; 829 } 830 831 /// \brief Starts the algorithm. 832 /// 833 /// Starts the algorithm. It runs augmenting phases until the 834 /// optimal solution reached. 835 /// 836 /// \param maxCardinality If the given value is true it will 837 /// calculate the maximum cardinality maximum matching instead of 838 /// the maximum matching. 839 void start(bool maxCardinality = false) { 840 while (augment(maxCardinality)) {} 841 } 842 843 /// \brief Runs the algorithm. 844 /// 845 /// It just initalize the algorithm and then start it. 846 /// 847 /// \param maxCardinality If the given value is true it will 848 /// calculate the maximum cardinality maximum matching instead of 849 /// the maximum matching. 850 void run(bool maxCardinality = false) { 851 init(); 852 start(maxCardinality); 853 } 854 855 /// @} 856 857 /// \name Query Functions 858 /// The result of the %Matching algorithm can be obtained using these 859 /// functions.\n 860 /// Before the use of these functions, 861 /// either run() or start() must be called. 862 863 ///@{ 864 865 /// \brief Gives back the potential in the NodeMap 866 /// 867 /// Gives back the potential in the NodeMap. The potential 868 /// is feasible if \f$ \pi(a) - \pi(b) - w(ab) = 0 \f$ for 869 /// each matching edges and \f$ \pi(a) - \pi(b) - w(ab) \ge 0 \f$ 870 /// for each edges. 871 template <typename PotentialMap> 872 void potential(PotentialMap& potential) { 873 for (ANodeIt it(*graph); it != INVALID; ++it) { 874 potential[it] = anode_potential[it]; 875 } 876 for (BNodeIt it(*graph); it != INVALID; ++it) { 877 potential[it] = bnode_potential[it]; 878 } 879 } 880 881 /// \brief Set true all matching uedge in the map. 882 /// 883 /// Set true all matching uedge in the map. It does not change the 884 /// value mapped to the other uedges. 885 /// \return The number of the matching edges. 886 template <typename MatchingMap> 887 int quickMatching(MatchingMap& matching) { 888 for (ANodeIt it(*graph); it != INVALID; ++it) { 889 if (anode_matching[it] != INVALID) { 890 matching[anode_matching[it]] = true; 891 } 892 } 893 return matching_size; 894 } 895 896 /// \brief Set true all matching uedge in the map and the others to false. 897 /// 898 /// Set true all matching uedge in the map and the others to false. 899 /// \return The number of the matching edges. 900 template <typename MatchingMap> 901 int matching(MatchingMap& matching) { 902 for (UEdgeIt it(*graph); it != INVALID; ++it) { 903 matching[it] = it == anode_matching[graph->aNode(it)]; 904 } 905 return matching_size; 906 } 907 908 909 /// \brief Return true if the given uedge is in the matching. 910 /// 911 /// It returns true if the given uedge is in the matching. 912 bool matchingEdge(const UEdge& edge) { 913 return anode_matching[graph->aNode(edge)] == edge; 914 } 915 916 /// \brief Returns the matching edge from the node. 917 /// 918 /// Returns the matching edge from the node. If there is not such 919 /// edge it gives back \c INVALID. 920 UEdge matchingEdge(const Node& node) { 921 if (graph->aNode(node)) { 922 return anode_matching[node]; 923 } else { 924 return bnode_matching[node]; 925 } 926 } 927 928 /// \brief Gives back the sum of weights of the matching edges. 929 /// 930 /// Gives back the sum of weights of the matching edges. 931 Value matchingValue() const { 932 return matching_value; 933 } 934 935 /// \brief Gives back the number of the matching edges. 936 /// 937 /// Gives back the number of the matching edges. 938 int matchingSize() const { 939 return matching_size; 940 } 941 942 /// @} 943 944 private: 945 946 void initStructures() { 947 if (!_heap_cross_ref) { 948 local_heap_cross_ref = true; 949 _heap_cross_ref = Traits::createHeapCrossRef(*graph); 950 } 951 if (!_heap) { 952 local_heap = true; 953 _heap = Traits::createHeap(*_heap_cross_ref); 954 } 955 } 956 957 void destroyStructures() { 958 if (local_heap_cross_ref) delete _heap_cross_ref; 959 if (local_heap) delete _heap; 960 } 961 962 963 private: 964 965 const BpUGraph *graph; 966 const WeightMap* weight; 967 968 ANodeMatchingMap anode_matching; 969 BNodeMatchingMap bnode_matching; 970 971 ANodePotentialMap anode_potential; 972 BNodePotentialMap bnode_potential; 973 974 Value matching_value; 975 int matching_size; 976 977 HeapCrossRef *_heap_cross_ref; 978 bool local_heap_cross_ref; 979 980 Heap *_heap; 981 bool local_heap; 982 983 }; 984 985 /// \brief Default traits class for minimum cost bipartite matching 986 /// algoritms. 987 /// 988 /// Default traits class for minimum cost bipartite matching 989 /// algoritms. 990 /// 991 /// \param _BpUGraph The bipartite undirected graph 992 /// type. 993 /// 994 /// \param _CostMap Type of cost map. 995 template <typename _BpUGraph, typename _CostMap> 996 struct MinCostMaxBipartiteMatchingDefaultTraits { 997 /// \brief The type of the cost of the undirected edges. 998 typedef typename _CostMap::Value Value; 999 1000 /// The undirected bipartite graph type the algorithm runs on. 1001 typedef _BpUGraph BpUGraph; 1002 1003 /// The map of the edges costs 1004 typedef _CostMap CostMap; 1005 1006 /// \brief The cross reference type used by heap. 1007 /// 1008 /// The cross reference type used by heap. 1009 /// Usually it is \c Graph::NodeMap<int>. 1010 typedef typename BpUGraph::template NodeMap<int> HeapCrossRef; 1011 1012 /// \brief Instantiates a HeapCrossRef. 1013 /// 1014 /// This function instantiates a \ref HeapCrossRef. 1015 /// \param graph is the graph, to which we would like to define the 1016 /// HeapCrossRef. 1017 static HeapCrossRef *createHeapCrossRef(const BpUGraph &graph) { 1018 return new HeapCrossRef(graph); 1019 } 1020 1021 /// \brief The heap type used by costed matching algorithms. 1022 /// 1023 /// The heap type used by costed matching algorithms. It should 1024 /// minimize the priorities and the heap's key type is the graph's 1025 /// anode graph's node. 1026 /// 1027 /// \sa BinHeap 1028 typedef BinHeap<typename BpUGraph::Node, Value, HeapCrossRef> Heap; 1029 1030 /// \brief Instantiates a Heap. 1031 /// 1032 /// This function instantiates a \ref Heap. 1033 /// \param crossref The cross reference of the heap. 1034 static Heap *createHeap(HeapCrossRef& crossref) { 1035 return new Heap(crossref); 1036 } 1037 1038 }; 1039 1040 1041 /// \ingroup matching 1042 /// 1043 /// \brief Bipartite Min Cost Matching algorithm 1044 /// 1045 /// This class implements the bipartite Min Cost Matching algorithm. 1046 /// It uses the successive shortest path algorithm to calculate the 1047 /// minimum cost maximum matching in the bipartite graph. The time 1048 /// complexity of the algorithm is \f$ O(ne\log(n)) \f$ with the 1049 /// default binary heap implementation but this can be improved to 1050 /// \f$ O(n^2\log(n)+ne) \f$ if we use fibonacci heaps. 1051 /// 1052 /// The algorithm also provides a potential function on the nodes 1053 /// which a dual solution of the matching algorithm and it can be 1054 /// used to proof the optimality of the given pimal solution. 1055 #ifdef DOXYGEN 1056 template <typename _BpUGraph, typename _CostMap, typename _Traits> 1057 #else 1058 template <typename _BpUGraph, 1059 typename _CostMap = typename _BpUGraph::template UEdgeMap<int>, 1060 typename _Traits = MinCostMaxBipartiteMatchingDefaultTraits<_BpUGraph, _CostMap> > 1061 #endif 1062 class MinCostMaxBipartiteMatching { 1063 public: 1064 1065 typedef _Traits Traits; 1066 typedef typename Traits::BpUGraph BpUGraph; 1067 typedef typename Traits::CostMap CostMap; 1068 typedef typename Traits::Value Value; 1069 1070 protected: 1071 1072 typedef typename Traits::HeapCrossRef HeapCrossRef; 1073 typedef typename Traits::Heap Heap; 1074 1075 1076 typedef typename BpUGraph::Node Node; 1077 typedef typename BpUGraph::ANodeIt ANodeIt; 1078 typedef typename BpUGraph::BNodeIt BNodeIt; 1079 typedef typename BpUGraph::UEdge UEdge; 1080 typedef typename BpUGraph::UEdgeIt UEdgeIt; 1081 typedef typename BpUGraph::IncEdgeIt IncEdgeIt; 1082 1083 typedef typename BpUGraph::template ANodeMap<UEdge> ANodeMatchingMap; 1084 typedef typename BpUGraph::template BNodeMap<UEdge> BNodeMatchingMap; 1085 1086 typedef typename BpUGraph::template ANodeMap<Value> ANodePotentialMap; 1087 typedef typename BpUGraph::template BNodeMap<Value> BNodePotentialMap; 1088 1089 1090 public: 1091 1092 /// \brief \ref Exception for uninitialized parameters. 1093 /// 1094 /// This error represents problems in the initialization 1095 /// of the parameters of the algorithms. 1096 class UninitializedParameter : public lemon::UninitializedParameter { 1097 public: 1098 virtual const char* exceptionName() const { 1099 return "lemon::MinCostMaxBipartiteMatching::UninitializedParameter"; 1100 } 1101 }; 1102 1103 ///\name Named template parameters 1104 1105 ///@{ 1106 1107 template <class H, class CR> 1108 struct DefHeapTraits : public Traits { 1109 typedef CR HeapCrossRef; 1110 typedef H Heap; 1111 static HeapCrossRef *createHeapCrossRef(const BpUGraph &) { 1112 throw UninitializedParameter(); 1113 } 1114 static Heap *createHeap(HeapCrossRef &) { 1115 throw UninitializedParameter(); 1116 } 1117 }; 1118 1119 /// \brief \ref named-templ-param "Named parameter" for setting heap 1120 /// and cross reference type 1121 /// 1122 /// \ref named-templ-param "Named parameter" for setting heap and cross 1123 /// reference type 1124 template <class H, class CR = typename BpUGraph::template NodeMap<int> > 1125 struct DefHeap 1126 : public MinCostMaxBipartiteMatching<BpUGraph, CostMap, 1127 DefHeapTraits<H, CR> > { 1128 typedef MinCostMaxBipartiteMatching<BpUGraph, CostMap, 1129 DefHeapTraits<H, CR> > Create; 1130 }; 1131 1132 template <class H, class CR> 1133 struct DefStandardHeapTraits : public Traits { 1134 typedef CR HeapCrossRef; 1135 typedef H Heap; 1136 static HeapCrossRef *createHeapCrossRef(const BpUGraph &graph) { 1137 return new HeapCrossRef(graph); 1138 } 1139 static Heap *createHeap(HeapCrossRef &crossref) { 1140 return new Heap(crossref); 1141 } 1142 }; 1143 1144 /// \brief \ref named-templ-param "Named parameter" for setting heap and 1145 /// cross reference type with automatic allocation 1146 /// 1147 /// \ref named-templ-param "Named parameter" for setting heap and cross 1148 /// reference type. It can allocate the heap and the cross reference 1149 /// object if the cross reference's constructor waits for the graph as 1150 /// parameter and the heap's constructor waits for the cross reference. 1151 template <class H, class CR = typename BpUGraph::template NodeMap<int> > 1152 struct DefStandardHeap 1153 : public MinCostMaxBipartiteMatching<BpUGraph, CostMap, 1154 DefStandardHeapTraits<H, CR> > { 1155 typedef MinCostMaxBipartiteMatching<BpUGraph, CostMap, 1156 DefStandardHeapTraits<H, CR> > 1157 Create; 1158 }; 1159 1160 ///@} 1161 1162 1163 /// \brief Constructor. 1164 /// 1165 /// Constructor of the algorithm. 1166 MinCostMaxBipartiteMatching(const BpUGraph& _graph, 1167 const CostMap& _cost) 1168 : graph(&_graph), cost(&_cost), 1169 anode_matching(_graph), bnode_matching(_graph), 1170 anode_potential(_graph), bnode_potential(_graph), 1171 _heap_cross_ref(0), local_heap_cross_ref(false), 1172 _heap(0), local_heap(0) {} 1173 1174 /// \brief Destructor. 1175 /// 1176 /// Destructor of the algorithm. 1177 ~MinCostMaxBipartiteMatching() { 1178 destroyStructures(); 1179 } 1180 1181 /// \brief Sets the heap and the cross reference used by algorithm. 1182 /// 1183 /// Sets the heap and the cross reference used by algorithm. 1184 /// If you don't use this function before calling \ref run(), 1185 /// it will allocate one. The destuctor deallocates this 1186 /// automatically allocated map, of course. 1187 /// \return \c (*this) 1188 MinCostMaxBipartiteMatching& heap(Heap& heap, HeapCrossRef &crossRef) { 1189 if(local_heap_cross_ref) { 1190 delete _heap_cross_ref; 1191 local_heap_cross_ref = false; 1192 } 1193 _heap_cross_ref = &crossRef; 1194 if(local_heap) { 1195 delete _heap; 1196 local_heap = false; 1197 } 1198 _heap = &heap; 1199 return *this; 1200 } 1201 1202 /// \name Execution control 1203 /// The simplest way to execute the algorithm is to use 1204 /// one of the member functions called \c run(). 1205 /// \n 1206 /// If you need more control on the execution, 1207 /// first you must call \ref init() or one alternative for it. 1208 /// Finally \ref start() will perform the matching computation or 1209 /// with step-by-step execution you can augment the solution. 1210 1211 /// @{ 1212 1213 /// \brief Initalize the data structures. 1214 /// 1215 /// It initalizes the data structures and creates an empty matching. 1216 void init() { 1217 initStructures(); 1218 for (ANodeIt it(*graph); it != INVALID; ++it) { 1219 anode_matching[it] = INVALID; 1220 anode_potential[it] = 0; 1221 } 1222 for (BNodeIt it(*graph); it != INVALID; ++it) { 1223 bnode_matching[it] = INVALID; 1224 bnode_potential[it] = 0; 1225 } 1226 matching_cost = 0; 1227 matching_size = 0; 1228 } 1229 1230 1231 /// \brief An augmenting phase of the costed matching algorithm 1232 /// 1233 /// It runs an augmenting phase of the matching algorithm. The 1234 /// phase finds the best augmenting path and augments only on this 1235 /// paths. 1236 /// 1237 /// The algorithm consists at most 1238 /// of \f$ O(n) \f$ phase and one phase is \f$ O(n\log(n)+e) \f$ 1239 /// long with Fibonacci heap or \f$ O((n+e)\log(n)) \f$ long 1240 /// with binary heap. 1241 bool augment() { 1242 1243 typename BpUGraph::template BNodeMap<Value> bdist(*graph); 1244 typename BpUGraph::template BNodeMap<UEdge> bpred(*graph, INVALID); 1245 1246 Node bestNode = INVALID; 1247 Value bestValue = 0; 1248 1249 _heap->clear(); 1250 for (ANodeIt it(*graph); it != INVALID; ++it) { 1251 (*_heap_cross_ref)[it] = Heap::PRE_HEAP; 1252 } 1253 1254 for (ANodeIt it(*graph); it != INVALID; ++it) { 1255 if (anode_matching[it] == INVALID) { 1256 _heap->push(it, 0); 1257 } 1258 } 1259 1260 while (!_heap->empty()) { 1261 Node anode = _heap->top(); 1262 Value avalue = _heap->prio(); 1263 _heap->pop(); 1264 for (IncEdgeIt jt(*graph, anode); jt != INVALID; ++jt) { 1265 if (jt == anode_matching[anode]) continue; 1266 Node bnode = graph->bNode(jt); 1267 Value bvalue = avalue + (*cost)[jt] + 1268 anode_potential[anode] - bnode_potential[bnode]; 1269 if (bpred[bnode] == INVALID || bvalue < bdist[bnode]) { 1270 bdist[bnode] = bvalue; 1271 bpred[bnode] = jt; 1272 } 1273 if (bnode_matching[bnode] != INVALID) { 1274 Node newanode = graph->aNode(bnode_matching[bnode]); 1275 switch (_heap->state(newanode)) { 1276 case Heap::PRE_HEAP: 1277 _heap->push(newanode, bvalue); 1278 break; 1279 case Heap::IN_HEAP: 1280 if (bvalue < (*_heap)[newanode]) { 1281 _heap->decrease(newanode, bvalue); 1282 } 1283 break; 1284 case Heap::POST_HEAP: 1285 break; 1286 } 1287 } else { 1288 if (bestNode == INVALID || 1289 bvalue + bnode_potential[bnode] < bestValue) { 1290 bestValue = bvalue + bnode_potential[bnode]; 1291 bestNode = bnode; 1292 } 1293 } 1294 } 1295 } 1296 1297 if (bestNode == INVALID) { 1298 return false; 1299 } 1300 1301 matching_cost += bestValue; 1302 ++matching_size; 1303 1304 for (BNodeIt it(*graph); it != INVALID; ++it) { 1305 if (bpred[it] != INVALID) { 1306 bnode_potential[it] += bdist[it]; 1307 } 1308 } 1309 for (ANodeIt it(*graph); it != INVALID; ++it) { 1310 if (anode_matching[it] != INVALID) { 1311 Node bnode = graph->bNode(anode_matching[it]); 1312 if (bpred[bnode] != INVALID) { 1313 anode_potential[it] += bdist[bnode]; 1314 } 1315 } 1316 } 1317 1318 while (bestNode != INVALID) { 1319 UEdge uedge = bpred[bestNode]; 1320 Node anode = graph->aNode(uedge); 1321 1322 bnode_matching[bestNode] = uedge; 1323 if (anode_matching[anode] != INVALID) { 1324 bestNode = graph->bNode(anode_matching[anode]); 1325 } else { 1326 bestNode = INVALID; 1327 } 1328 anode_matching[anode] = uedge; 1329 } 1330 1331 1332 return true; 1333 } 1334 1335 /// \brief Starts the algorithm. 1336 /// 1337 /// Starts the algorithm. It runs augmenting phases until the 1338 /// optimal solution reached. 1339 void start() { 1340 while (augment()) {} 1341 } 1342 1343 /// \brief Runs the algorithm. 1344 /// 1345 /// It just initalize the algorithm and then start it. 1346 void run() { 1347 init(); 1348 start(); 1349 } 1350 1351 /// @} 1352 1353 /// \name Query Functions 1354 /// The result of the %Matching algorithm can be obtained using these 1355 /// functions.\n 1356 /// Before the use of these functions, 1357 /// either run() or start() must be called. 1358 1359 ///@{ 1360 1361 /// \brief Gives back the potential in the NodeMap 1362 /// 1363 /// Gives back the potential in the NodeMap. The potential 1364 /// is feasible if \f$ \pi(a) - \pi(b) + w(ab) = 0 \f$ for 1365 /// each matching edges and \f$ \pi(a) - \pi(b) + w(ab) \ge 0 \f$ 1366 /// for each edges. 1367 template <typename PotentialMap> 1368 void potential(PotentialMap& potential) { 1369 for (ANodeIt it(*graph); it != INVALID; ++it) { 1370 potential[it] = anode_potential[it]; 1371 } 1372 for (BNodeIt it(*graph); it != INVALID; ++it) { 1373 potential[it] = bnode_potential[it]; 1374 } 1375 } 1376 1377 /// \brief Set true all matching uedge in the map. 1378 /// 1379 /// Set true all matching uedge in the map. It does not change the 1380 /// value mapped to the other uedges. 1381 /// \return The number of the matching edges. 1382 template <typename MatchingMap> 1383 int quickMatching(MatchingMap& matching) { 1384 for (ANodeIt it(*graph); it != INVALID; ++it) { 1385 if (anode_matching[it] != INVALID) { 1386 matching[anode_matching[it]] = true; 1387 } 1388 } 1389 return matching_size; 1390 } 1391 1392 /// \brief Set true all matching uedge in the map and the others to false. 1393 /// 1394 /// Set true all matching uedge in the map and the others to false. 1395 /// \return The number of the matching edges. 1396 template <typename MatchingMap> 1397 int matching(MatchingMap& matching) { 1398 for (UEdgeIt it(*graph); it != INVALID; ++it) { 1399 matching[it] = it == anode_matching[graph->aNode(it)]; 1400 } 1401 return matching_size; 1402 } 1403 1404 1405 /// \brief Return true if the given uedge is in the matching. 1406 /// 1407 /// It returns true if the given uedge is in the matching. 1408 bool matchingEdge(const UEdge& edge) { 1409 return anode_matching[graph->aNode(edge)] == edge; 1410 } 1411 1412 /// \brief Returns the matching edge from the node. 1413 /// 1414 /// Returns the matching edge from the node. If there is not such 1415 /// edge it gives back \c INVALID. 1416 UEdge matchingEdge(const Node& node) { 1417 if (graph->aNode(node)) { 1418 return anode_matching[node]; 1419 } else { 1420 return bnode_matching[node]; 1421 } 1422 } 1423 1424 /// \brief Gives back the sum of costs of the matching edges. 1425 /// 1426 /// Gives back the sum of costs of the matching edges. 1427 Value matchingCost() const { 1428 return matching_cost; 1429 } 1430 1431 /// \brief Gives back the number of the matching edges. 1432 /// 1433 /// Gives back the number of the matching edges. 1434 int matchingSize() const { 1435 return matching_size; 1436 } 1437 1438 /// @} 1439 1440 private: 1441 1442 void initStructures() { 1443 if (!_heap_cross_ref) { 1444 local_heap_cross_ref = true; 1445 _heap_cross_ref = Traits::createHeapCrossRef(*graph); 1446 } 1447 if (!_heap) { 1448 local_heap = true; 1449 _heap = Traits::createHeap(*_heap_cross_ref); 1450 } 1451 } 1452 1453 void destroyStructures() { 1454 if (local_heap_cross_ref) delete _heap_cross_ref; 1455 if (local_heap) delete _heap; 1456 } 1457 1458 1459 private: 1460 1461 const BpUGraph *graph; 1462 const CostMap* cost; 1463 1464 ANodeMatchingMap anode_matching; 1465 BNodeMatchingMap bnode_matching; 1466 1467 ANodePotentialMap anode_potential; 1468 BNodePotentialMap bnode_potential; 1469 1470 Value matching_cost; 1471 int matching_size; 1472 1473 HeapCrossRef *_heap_cross_ref; 1474 bool local_heap_cross_ref; 1475 1476 Heap *_heap; 1477 bool local_heap; 1478 1479 }; 1480 464 1481 } 465 1482 -
test/bipartite_matching_test.cc
r2040 r2051 9 9 10 10 #include <lemon/graph_utils.h> 11 12 #include <lemon/maps.h> 11 13 12 14 #include "test_tools.h" … … 25 27 int m = argc > 2 ? atoi(argv[2]) : 100; 26 28 int e = argc > 3 ? atoi(argv[3]) : (int)((n+m) * log(n+m)); 29 int c = argc > 4 ? atoi(argv[4]) : 100; 30 31 Graph::UEdgeMap<int> weight(graph); 32 33 int max_cardinality; 34 int max_weight; 35 int max_cardinality_max_weight; 27 36 28 37 for (int i = 0; i < n; ++i) { … … 37 46 Node aNode = aNodes[urandom(n)]; 38 47 Node bNode = bNodes[urandom(m)]; 39 graph.addEdge(aNode, bNode); 48 UEdge uedge = graph.addEdge(aNode, bNode); 49 weight[uedge] = urandom(c); 40 50 } 41 51 … … 58 68 for (IncEdgeIt jt(graph, it); jt != INVALID; ++jt) { 59 69 if (mm[jt]) ++num; 60 } 61 check(num <= 1, "INVALID PRIMAL"); 62 } 70 } 71 check(num <= 1, "INVALID PRIMAL"); 72 } 73 max_cardinality = bpmatch.matchingSize(); 63 74 } 64 75 … … 112 123 113 124 { 114 SwapBpUGraphAdaptor<Graph> swapgraph(graph); 115 MaxBipartiteMatching<SwapBpUGraphAdaptor<Graph> > bpmatch(swapgraph); 125 MaxWeightedBipartiteMatching<Graph> bpmatch(graph, weight); 126 127 bpmatch.init(); 128 129 max_weight = 0; 130 while (bpmatch.augment(true)) { 131 132 Graph::UEdgeMap<bool> mm(graph); 133 Graph::NodeMap<int> pm(graph); 134 135 bpmatch.matching(mm); 136 bpmatch.potential(pm); 137 138 for (UEdgeIt it(graph); it != INVALID; ++it) { 139 if (mm[it]) { 140 check(pm[graph.aNode(it)] - weight[it] - pm[graph.bNode(it)] == 0, 141 "INVALID DUAL"); 142 } else { 143 check(pm[graph.aNode(it)] - weight[it] - pm[graph.bNode(it)] >= 0, 144 "INVALID DUAL"); 145 } 146 } 147 148 for (ANodeIt it(graph); it != INVALID; ++it) { 149 int num = 0; 150 for (IncEdgeIt jt(graph, it); jt != INVALID; ++jt) { 151 if (mm[jt]) ++num; 152 } 153 check(num <= 1, "INVALID PRIMAL"); 154 } 155 if (bpmatch.matchingValue() > max_weight) { 156 max_weight = bpmatch.matchingValue(); 157 } 158 } 159 } 160 161 { 162 MaxWeightedBipartiteMatching<Graph> bpmatch(graph, weight); 116 163 117 164 bpmatch.run(); 118 119 Graph::UEdgeMap<bool> mm(graph); 120 Graph::NodeMap<bool> cs(graph); 121 122 check(bpmatch.coverSet(cs) == bpmatch.matching(mm), "INVALID PRIMAL-DUAL"); 123 124 for (UEdgeIt it(graph); it != INVALID; ++it) { 125 check(cs[graph.aNode(it)] || cs[graph.bNode(it)], "INVALID DUAL"); 126 } 127 128 for (ANodeIt it(graph); it != INVALID; ++it) { 129 int num = 0; 130 for (IncEdgeIt jt(graph, it); jt != INVALID; ++jt) { 131 if (mm[jt]) ++num; 132 } 133 check(num <= 1, "INVALID PRIMAL"); 134 } 165 166 Graph::UEdgeMap<bool> mm(graph); 167 Graph::NodeMap<int> pm(graph); 168 169 bpmatch.matching(mm); 170 bpmatch.potential(pm); 171 172 for (UEdgeIt it(graph); it != INVALID; ++it) { 173 if (mm[it]) { 174 check(pm[graph.aNode(it)] - weight[it] - pm[graph.bNode(it)] == 0, 175 "INVALID DUAL"); 176 } else { 177 check(pm[graph.aNode(it)] - weight[it] - pm[graph.bNode(it)] >= 0, 178 "INVALID DUAL"); 179 } 180 } 181 182 for (ANodeIt it(graph); it != INVALID; ++it) { 183 int num = 0; 184 for (IncEdgeIt jt(graph, it); jt != INVALID; ++jt) { 185 if (mm[jt]) ++num; 186 } 187 check(num <= 1, "INVALID PRIMAL"); 188 } 189 check(max_weight == bpmatch.matchingValue(), "WRONG WEIGHT"); 190 } 191 192 { 193 MaxWeightedBipartiteMatching<Graph> bpmatch(graph, weight); 194 195 bpmatch.run(true); 196 197 Graph::UEdgeMap<bool> mm(graph); 198 Graph::NodeMap<int> pm(graph); 199 200 bpmatch.matching(mm); 201 bpmatch.potential(pm); 202 203 for (UEdgeIt it(graph); it != INVALID; ++it) { 204 if (mm[it]) { 205 check(pm[graph.aNode(it)] - weight[it] - pm[graph.bNode(it)] == 0, 206 "INVALID DUAL"); 207 } else { 208 check(pm[graph.aNode(it)] - weight[it] - pm[graph.bNode(it)] >= 0, 209 "INVALID DUAL"); 210 } 211 } 212 213 for (ANodeIt it(graph); it != INVALID; ++it) { 214 int num = 0; 215 for (IncEdgeIt jt(graph, it); jt != INVALID; ++jt) { 216 if (mm[jt]) ++num; 217 } 218 check(num <= 1, "INVALID PRIMAL"); 219 } 220 check(max_cardinality == bpmatch.matchingSize(), "WRONG SIZE"); 221 max_cardinality_max_weight = bpmatch.matchingValue(); 222 223 } 224 225 { 226 Graph::UEdgeMap<int> cost(graph); 227 228 cost = subMap(constMap<UEdge>(c), weight); 229 230 MinCostMaxBipartiteMatching<Graph> bpmatch(graph, cost); 231 232 bpmatch.run(); 233 234 Graph::UEdgeMap<bool> mm(graph); 235 Graph::NodeMap<int> pm(graph); 236 237 bpmatch.matching(mm); 238 bpmatch.potential(pm); 239 240 for (UEdgeIt it(graph); it != INVALID; ++it) { 241 if (mm[it]) { 242 check(pm[graph.aNode(it)] + cost[it] - pm[graph.bNode(it)] == 0, 243 "INVALID DUAL"); 244 } else { 245 check(pm[graph.aNode(it)] + cost[it] - pm[graph.bNode(it)] >= 0, 246 "INVALID DUAL"); 247 } 248 } 249 250 for (ANodeIt it(graph); it != INVALID; ++it) { 251 int num = 0; 252 for (IncEdgeIt jt(graph, it); jt != INVALID; ++jt) { 253 if (mm[jt]) ++num; 254 } 255 check(num <= 1, "INVALID PRIMAL"); 256 } 257 258 check(max_cardinality == bpmatch.matchingSize(), "WRONG SIZE"); 259 check(max_cardinality * c - max_cardinality_max_weight 260 == bpmatch.matchingCost(), "WRONG SIZE"); 261 135 262 } 136 263
Note: See TracChangeset
for help on using the changeset viewer.