... | ... |
@@ -17,6 +17,6 @@ |
17 | 17 |
*/ |
18 | 18 |
|
19 |
#ifndef LEMON_MAX_MATCHING_H |
|
20 |
#define LEMON_MAX_MATCHING_H |
|
19 |
#ifndef LEMON_MATCHING_H |
|
20 |
#define LEMON_MATCHING_H |
|
21 | 21 |
|
22 | 22 |
#include <vector> |
... | ... |
@@ -42,5 +42,5 @@ |
42 | 42 |
/// This class implements Edmonds' alternating forest matching algorithm |
43 | 43 |
/// for finding a maximum cardinality matching in a general undirected graph. |
44 |
/// It can be started from an arbitrary initial matching |
|
44 |
/// It can be started from an arbitrary initial matching |
|
45 | 45 |
/// (the default is the empty one). |
46 | 46 |
/// |
... | ... |
@@ -70,9 +70,9 @@ |
70 | 70 |
///\brief Status constants for Gallai-Edmonds decomposition. |
71 | 71 |
/// |
72 |
///These constants are used for indicating the Gallai-Edmonds |
|
72 |
///These constants are used for indicating the Gallai-Edmonds |
|
73 | 73 |
///decomposition of a graph. The nodes with status \c EVEN (or \c D) |
74 | 74 |
///induce a subgraph with factor-critical components, the nodes with |
75 | 75 |
///status \c ODD (or \c A) form the canonical barrier, and the nodes |
76 |
///with status \c MATCHED (or \c C) induce a subgraph having a |
|
76 |
///with status \c MATCHED (or \c C) induce a subgraph having a |
|
77 | 77 |
///perfect matching. |
78 | 78 |
enum Status { |
... | ... |
@@ -513,5 +513,5 @@ |
513 | 513 |
} |
514 | 514 |
|
515 |
/// \brief Start Edmonds' algorithm with a heuristic improvement |
|
515 |
/// \brief Start Edmonds' algorithm with a heuristic improvement |
|
516 | 516 |
/// for dense graphs |
517 | 517 |
/// |
... | ... |
@@ -535,6 +535,6 @@ |
535 | 535 |
/// \brief Run Edmonds' algorithm |
536 | 536 |
/// |
537 |
/// This function runs Edmonds' algorithm. An additional heuristic of |
|
538 |
/// postponing shrinks is used for relatively dense graphs |
|
537 |
/// This function runs Edmonds' algorithm. An additional heuristic of |
|
538 |
/// postponing shrinks is used for relatively dense graphs |
|
539 | 539 |
/// (for which <tt>m>=2*n</tt> holds). |
540 | 540 |
void run() { |
... | ... |
@@ -557,5 +557,5 @@ |
557 | 557 |
/// \brief Return the size (cardinality) of the matching. |
558 | 558 |
/// |
559 |
/// This function returns the size (cardinality) of the current matching. |
|
559 |
/// This function returns the size (cardinality) of the current matching. |
|
560 | 560 |
/// After run() it returns the size of the maximum matching in the graph. |
561 | 561 |
int matchingSize() const { |
... | ... |
@@ -571,5 +571,5 @@ |
571 | 571 |
/// \brief Return \c true if the given edge is in the matching. |
572 | 572 |
/// |
573 |
/// This function returns \c true if the given edge is in the current |
|
573 |
/// This function returns \c true if the given edge is in the current |
|
574 | 574 |
/// matching. |
575 | 575 |
bool matching(const Edge& edge) const { |
... | ... |
@@ -580,5 +580,5 @@ |
580 | 580 |
/// |
581 | 581 |
/// This function returns the matching arc (or edge) incident to the |
582 |
/// given node in the current matching or \c INVALID if the node is |
|
582 |
/// given node in the current matching or \c INVALID if the node is |
|
583 | 583 |
/// not covered by the matching. |
584 | 584 |
Arc matching(const Node& n) const { |
... | ... |
@@ -596,5 +596,5 @@ |
596 | 596 |
/// \brief Return the mate of the given node. |
597 | 597 |
/// |
598 |
/// This function returns the mate of the given node in the current |
|
598 |
/// This function returns the mate of the given node in the current |
|
599 | 599 |
/// matching or \c INVALID if the node is not covered by the matching. |
600 | 600 |
Node mate(const Node& n) const { |
... | ... |
@@ -606,5 +606,5 @@ |
606 | 606 |
|
607 | 607 |
/// \name Dual Solution |
608 |
/// Functions to get the dual solution, i.e. the Gallai-Edmonds |
|
608 |
/// Functions to get the dual solution, i.e. the Gallai-Edmonds |
|
609 | 609 |
/// decomposition. |
610 | 610 |
|
... | ... |
@@ -649,6 +649,6 @@ |
649 | 649 |
/// \f$O(nm\log n)\f$ time complexity. |
650 | 650 |
/// |
651 |
/// The maximum weighted matching problem is to find a subset of the |
|
652 |
/// edges in an undirected graph with maximum overall weight for which |
|
651 |
/// The maximum weighted matching problem is to find a subset of the |
|
652 |
/// edges in an undirected graph with maximum overall weight for which |
|
653 | 653 |
/// each node has at most one incident edge. |
654 | 654 |
/// It can be formulated with the following linear program. |
... | ... |
@@ -674,14 +674,14 @@ |
674 | 674 |
\frac{\vert B \vert - 1}{2}z_B\f] */ |
675 | 675 |
/// |
676 |
/// The algorithm can be executed with the run() function. |
|
676 |
/// The algorithm can be executed with the run() function. |
|
677 | 677 |
/// After it the matching (the primal solution) and the dual solution |
678 |
/// can be obtained using the query functions and the |
|
679 |
/// \ref MaxWeightedMatching::BlossomIt "BlossomIt" nested class, |
|
680 |
/// |
|
678 |
/// can be obtained using the query functions and the |
|
679 |
/// \ref MaxWeightedMatching::BlossomIt "BlossomIt" nested class, |
|
680 |
/// which is able to iterate on the nodes of a blossom. |
|
681 | 681 |
/// If the value type is integer, then the dual solution is multiplied |
682 | 682 |
/// by \ref MaxWeightedMatching::dualScale "4". |
683 | 683 |
/// |
684 | 684 |
/// \tparam GR The undirected graph type the algorithm runs on. |
685 |
/// \tparam WM The type edge weight map. The default type is |
|
685 |
/// \tparam WM The type edge weight map. The default type is |
|
686 | 686 |
/// \ref concepts::Graph::EdgeMap "GR::EdgeMap<int>". |
687 | 687 |
#ifdef DOXYGEN |
... | ... |
@@ -746,5 +746,5 @@ |
746 | 746 |
|
747 | 747 |
enum Status { |
748 |
EVEN = -1, MATCHED = 0, ODD = 1 |
|
748 |
EVEN = -1, MATCHED = 0, ODD = 1 |
|
749 | 749 |
}; |
750 | 750 |
|
... | ... |
@@ -845,7 +845,4 @@ |
845 | 845 |
|
846 | 846 |
void destroyStructures() { |
847 |
_node_num = countNodes(_graph); |
|
848 |
_blossom_num = _node_num * 3 / 2; |
|
849 |
|
|
850 | 847 |
if (_matching) { |
851 | 848 |
delete _matching; |
... | ... |
@@ -923,8 +920,4 @@ |
923 | 920 |
_delta3->push(e, rw / 2); |
924 | 921 |
} |
925 |
} else if ((*_blossom_data)[vb].status == UNMATCHED) { |
|
926 |
if (_delta3->state(e) != _delta3->IN_HEAP) { |
|
927 |
_delta3->push(e, rw); |
|
928 |
} |
|
929 | 922 |
} else { |
930 | 923 |
typename std::map<int, Arc>::iterator it = |
... | ... |
@@ -950,200 +943,4 @@ |
950 | 943 |
(*_blossom_data)[vb].offset); |
951 | 944 |
} else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) - |
952 |
(*_blossom_data)[vb].offset){ |
|
953 |
_delta2->decrease(vb, _blossom_set->classPrio(vb) - |
|
954 |
(*_blossom_data)[vb].offset); |
|
955 |
} |
|
956 |
} |
|
957 |
} |
|
958 |
} |
|
959 |
} |
|
960 |
} |
|
961 |
(*_blossom_data)[blossom].offset = 0; |
|
962 |
} |
|
963 |
|
|
964 |
void matchedToOdd(int blossom) { |
|
965 |
if (_delta2->state(blossom) == _delta2->IN_HEAP) { |
|
966 |
_delta2->erase(blossom); |
|
967 |
} |
|
968 |
(*_blossom_data)[blossom].offset += _delta_sum; |
|
969 |
if (!_blossom_set->trivial(blossom)) { |
|
970 |
_delta4->push(blossom, (*_blossom_data)[blossom].pot / 2 + |
|
971 |
(*_blossom_data)[blossom].offset); |
|
972 |
} |
|
973 |
} |
|
974 |
|
|
975 |
void evenToMatched(int blossom, int tree) { |
|
976 |
if (!_blossom_set->trivial(blossom)) { |
|
977 |
(*_blossom_data)[blossom].pot += 2 * _delta_sum; |
|
978 |
} |
|
979 |
|
|
980 |
for (typename BlossomSet::ItemIt n(*_blossom_set, blossom); |
|
981 |
n != INVALID; ++n) { |
|
982 |
int ni = (*_node_index)[n]; |
|
983 |
(*_node_data)[ni].pot -= _delta_sum; |
|
984 |
|
|
985 |
_delta1->erase(n); |
|
986 |
|
|
987 |
for (InArcIt e(_graph, n); e != INVALID; ++e) { |
|
988 |
Node v = _graph.source(e); |
|
989 |
int vb = _blossom_set->find(v); |
|
990 |
int vi = (*_node_index)[v]; |
|
991 |
|
|
992 |
Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot - |
|
993 |
dualScale * _weight[e]; |
|
994 |
|
|
995 |
if (vb == blossom) { |
|
996 |
if (_delta3->state(e) == _delta3->IN_HEAP) { |
|
997 |
_delta3->erase(e); |
|
998 |
} |
|
999 |
} else if ((*_blossom_data)[vb].status == EVEN) { |
|
1000 |
|
|
1001 |
if (_delta3->state(e) == _delta3->IN_HEAP) { |
|
1002 |
_delta3->erase(e); |
|
1003 |
} |
|
1004 |
|
|
1005 |
int vt = _tree_set->find(vb); |
|
1006 |
|
|
1007 |
if (vt != tree) { |
|
1008 |
|
|
1009 |
Arc r = _graph.oppositeArc(e); |
|
1010 |
|
|
1011 |
typename std::map<int, Arc>::iterator it = |
|
1012 |
(*_node_data)[ni].heap_index.find(vt); |
|
1013 |
|
|
1014 |
if (it != (*_node_data)[ni].heap_index.end()) { |
|
1015 |
if ((*_node_data)[ni].heap[it->second] > rw) { |
|
1016 |
(*_node_data)[ni].heap.replace(it->second, r); |
|
1017 |
(*_node_data)[ni].heap.decrease(r, rw); |
|
1018 |
it->second = r; |
|
1019 |
} |
|
1020 |
} else { |
|
1021 |
(*_node_data)[ni].heap.push(r, rw); |
|
1022 |
(*_node_data)[ni].heap_index.insert(std::make_pair(vt, r)); |
|
1023 |
} |
|
1024 |
|
|
1025 |
if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) { |
|
1026 |
_blossom_set->decrease(n, (*_node_data)[ni].heap.prio()); |
|
1027 |
|
|
1028 |
if (_delta2->state(blossom) != _delta2->IN_HEAP) { |
|
1029 |
_delta2->push(blossom, _blossom_set->classPrio(blossom) - |
|
1030 |
(*_blossom_data)[blossom].offset); |
|
1031 |
} else if ((*_delta2)[blossom] > |
|
1032 |
_blossom_set->classPrio(blossom) - |
|
1033 |
(*_blossom_data)[blossom].offset){ |
|
1034 |
_delta2->decrease(blossom, _blossom_set->classPrio(blossom) - |
|
1035 |
(*_blossom_data)[blossom].offset); |
|
1036 |
} |
|
1037 |
} |
|
1038 |
} |
|
1039 |
|
|
1040 |
} else if ((*_blossom_data)[vb].status == UNMATCHED) { |
|
1041 |
if (_delta3->state(e) == _delta3->IN_HEAP) { |
|
1042 |
_delta3->erase(e); |
|
1043 |
} |
|
1044 |
} else { |
|
1045 |
|
|
1046 |
typename std::map<int, Arc>::iterator it = |
|
1047 |
(*_node_data)[vi].heap_index.find(tree); |
|
1048 |
|
|
1049 |
if (it != (*_node_data)[vi].heap_index.end()) { |
|
1050 |
(*_node_data)[vi].heap.erase(it->second); |
|
1051 |
(*_node_data)[vi].heap_index.erase(it); |
|
1052 |
if ((*_node_data)[vi].heap.empty()) { |
|
1053 |
_blossom_set->increase(v, std::numeric_limits<Value>::max()); |
|
1054 |
} else if ((*_blossom_set)[v] < (*_node_data)[vi].heap.prio()) { |
|
1055 |
_blossom_set->increase(v, (*_node_data)[vi].heap.prio()); |
|
1056 |
} |
|
1057 |
|
|
1058 |
if ((*_blossom_data)[vb].status == MATCHED) { |
|
1059 |
if (_blossom_set->classPrio(vb) == |
|
1060 |
std::numeric_limits<Value>::max()) { |
|
1061 |
_delta2->erase(vb); |
|
1062 |
} else if ((*_delta2)[vb] < _blossom_set->classPrio(vb) - |
|
1063 |
(*_blossom_data)[vb].offset) { |
|
1064 |
_delta2->increase(vb, _blossom_set->classPrio(vb) - |
|
1065 |
(*_blossom_data)[vb].offset); |
|
1066 |
} |
|
1067 |
} |
|
1068 |
} |
|
1069 |
} |
|
1070 |
} |
|
1071 |
} |
|
1072 |
} |
|
1073 |
|
|
1074 |
void oddToMatched(int blossom) { |
|
1075 |
(*_blossom_data)[blossom].offset -= _delta_sum; |
|
1076 |
|
|
1077 |
if (_blossom_set->classPrio(blossom) != |
|
1078 |
std::numeric_limits<Value>::max()) { |
|
1079 |
_delta2->push(blossom, _blossom_set->classPrio(blossom) - |
|
1080 |
(*_blossom_data)[blossom].offset); |
|
1081 |
} |
|
1082 |
|
|
1083 |
if (!_blossom_set->trivial(blossom)) { |
|
1084 |
_delta4->erase(blossom); |
|
1085 |
} |
|
1086 |
} |
|
1087 |
|
|
1088 |
void oddToEven(int blossom, int tree) { |
|
1089 |
if (!_blossom_set->trivial(blossom)) { |
|
1090 |
_delta4->erase(blossom); |
|
1091 |
(*_blossom_data)[blossom].pot -= |
|
1092 |
2 * (2 * _delta_sum - (*_blossom_data)[blossom].offset); |
|
1093 |
} |
|
1094 |
|
|
1095 |
for (typename BlossomSet::ItemIt n(*_blossom_set, blossom); |
|
1096 |
n != INVALID; ++n) { |
|
1097 |
int ni = (*_node_index)[n]; |
|
1098 |
|
|
1099 |
_blossom_set->increase(n, std::numeric_limits<Value>::max()); |
|
1100 |
|
|
1101 |
(*_node_data)[ni].heap.clear(); |
|
1102 |
(*_node_data)[ni].heap_index.clear(); |
|
1103 |
(*_node_data)[ni].pot += |
|
1104 |
2 * _delta_sum - (*_blossom_data)[blossom].offset; |
|
1105 |
|
|
1106 |
_delta1->push(n, (*_node_data)[ni].pot); |
|
1107 |
|
|
1108 |
for (InArcIt e(_graph, n); e != INVALID; ++e) { |
|
1109 |
Node v = _graph.source(e); |
|
1110 |
int vb = _blossom_set->find(v); |
|
1111 |
int vi = (*_node_index)[v]; |
|
1112 |
|
|
1113 |
Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot - |
|
1114 |
dualScale * _weight[e]; |
|
1115 |
|
|
1116 |
if ((*_blossom_data)[vb].status == EVEN) { |
|
1117 |
if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) { |
|
1118 |
_delta3->push(e, rw / 2); |
|
1119 |
} |
|
1120 |
} else if ((*_blossom_data)[vb].status == UNMATCHED) { |
|
1121 |
if (_delta3->state(e) != _delta3->IN_HEAP) { |
|
1122 |
_delta3->push(e, rw); |
|
1123 |
} |
|
1124 |
} else { |
|
1125 |
|
|
1126 |
typename std::map<int, Arc>::iterator it = |
|
1127 |
(*_node_data)[vi].heap_index.find(tree); |
|
1128 |
|
|
1129 |
if (it != (*_node_data)[vi].heap_index.end()) { |
|
1130 |
if ((*_node_data)[vi].heap[it->second] > rw) { |
|
1131 |
(*_node_data)[vi].heap.replace(it->second, e); |
|
1132 |
(*_node_data)[vi].heap.decrease(e, rw); |
|
1133 |
it->second = e; |
|
1134 |
} |
|
1135 |
} else { |
|
1136 |
(*_node_data)[vi].heap.push(e, rw); |
|
1137 |
(*_node_data)[vi].heap_index.insert(std::make_pair(tree, e)); |
|
1138 |
} |
|
1139 |
|
|
1140 |
if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) { |
|
1141 |
_blossom_set->decrease(v, (*_node_data)[vi].heap.prio()); |
|
1142 |
|
|
1143 |
if ((*_blossom_data)[vb].status == MATCHED) { |
|
1144 |
if (_delta2->state(vb) != _delta2->IN_HEAP) { |
|
1145 |
_delta2->push(vb, _blossom_set->classPrio(vb) - |
|
1146 |
(*_blossom_data)[vb].offset); |
|
1147 |
} else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) - |
|
1148 | 945 |
(*_blossom_data)[vb].offset) { |
1149 | 946 |
_delta2->decrease(vb, _blossom_set->classPrio(vb) - |
... | ... |
@@ -1158,21 +955,29 @@ |
1158 | 955 |
} |
1159 | 956 |
|
1160 |
|
|
1161 |
void matchedToUnmatched(int blossom) { |
|
957 |
void matchedToOdd(int blossom) { |
|
1162 | 958 |
if (_delta2->state(blossom) == _delta2->IN_HEAP) { |
1163 | 959 |
_delta2->erase(blossom); |
1164 | 960 |
} |
961 |
(*_blossom_data)[blossom].offset += _delta_sum; |
|
962 |
if (!_blossom_set->trivial(blossom)) { |
|
963 |
_delta4->push(blossom, (*_blossom_data)[blossom].pot / 2 + |
|
964 |
(*_blossom_data)[blossom].offset); |
|
965 |
} |
|
966 |
} |
|
967 |
|
|
968 |
void evenToMatched(int blossom, int tree) { |
|
969 |
if (!_blossom_set->trivial(blossom)) { |
|
970 |
(*_blossom_data)[blossom].pot += 2 * _delta_sum; |
|
971 |
} |
|
1165 | 972 |
|
1166 | 973 |
for (typename BlossomSet::ItemIt n(*_blossom_set, blossom); |
1167 | 974 |
n != INVALID; ++n) { |
1168 | 975 |
int ni = (*_node_index)[n]; |
1169 |
|
|
1170 |
_blossom_set->increase(n, std::numeric_limits<Value>::max()); |
|
1171 |
|
|
1172 |
(*_node_data)[ni].heap.clear(); |
|
1173 |
(*_node_data)[ni].heap_index.clear(); |
|
1174 |
|
|
1175 |
for (OutArcIt e(_graph, n); e != INVALID; ++e) { |
|
1176 |
Node v = _graph.target(e); |
|
976 |
(*_node_data)[ni].pot -= _delta_sum; |
|
977 |
|
|
978 |
_delta1->erase(n); |
|
979 |
|
|
980 |
for (InArcIt e(_graph, n); e != INVALID; ++e) { |
|
981 |
Node v = _graph.source(e); |
|
1177 | 982 |
int vb = _blossom_set->find(v); |
1178 | 983 |
int vi = (*_node_index)[v]; |
... | ... |
@@ -1181,7 +986,72 @@ |
1181 | 986 |
dualScale * _weight[e]; |
1182 | 987 |
|
1183 |
if ((*_blossom_data)[vb].status == EVEN) { |
|
1184 |
if (_delta3->state(e) != _delta3->IN_HEAP) { |
|
1185 |
|
|
988 |
if (vb == blossom) { |
|
989 |
if (_delta3->state(e) == _delta3->IN_HEAP) { |
|
990 |
_delta3->erase(e); |
|
991 |
} |
|
992 |
} else if ((*_blossom_data)[vb].status == EVEN) { |
|
993 |
|
|
994 |
if (_delta3->state(e) == _delta3->IN_HEAP) { |
|
995 |
_delta3->erase(e); |
|
996 |
} |
|
997 |
|
|
998 |
int vt = _tree_set->find(vb); |
|
999 |
|
|
1000 |
if (vt != tree) { |
|
1001 |
|
|
1002 |
Arc r = _graph.oppositeArc(e); |
|
1003 |
|
|
1004 |
typename std::map<int, Arc>::iterator it = |
|
1005 |
(*_node_data)[ni].heap_index.find(vt); |
|
1006 |
|
|
1007 |
if (it != (*_node_data)[ni].heap_index.end()) { |
|
1008 |
if ((*_node_data)[ni].heap[it->second] > rw) { |
|
1009 |
(*_node_data)[ni].heap.replace(it->second, r); |
|
1010 |
(*_node_data)[ni].heap.decrease(r, rw); |
|
1011 |
it->second = r; |
|
1012 |
} |
|
1013 |
} else { |
|
1014 |
(*_node_data)[ni].heap.push(r, rw); |
|
1015 |
(*_node_data)[ni].heap_index.insert(std::make_pair(vt, r)); |
|
1016 |
} |
|
1017 |
|
|
1018 |
if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) { |
|
1019 |
_blossom_set->decrease(n, (*_node_data)[ni].heap.prio()); |
|
1020 |
|
|
1021 |
if (_delta2->state(blossom) != _delta2->IN_HEAP) { |
|
1022 |
_delta2->push(blossom, _blossom_set->classPrio(blossom) - |
|
1023 |
(*_blossom_data)[blossom].offset); |
|
1024 |
} else if ((*_delta2)[blossom] > |
|
1025 |
_blossom_set->classPrio(blossom) - |
|
1026 |
(*_blossom_data)[blossom].offset){ |
|
1027 |
_delta2->decrease(blossom, _blossom_set->classPrio(blossom) - |
|
1028 |
(*_blossom_data)[blossom].offset); |
|
1029 |
} |
|
1030 |
} |
|
1031 |
} |
|
1032 |
} else { |
|
1033 |
|
|
1034 |
typename std::map<int, Arc>::iterator it = |
|
1035 |
(*_node_data)[vi].heap_index.find(tree); |
|
1036 |
|
|
1037 |
if (it != (*_node_data)[vi].heap_index.end()) { |
|
1038 |
(*_node_data)[vi].heap.erase(it->second); |
|
1039 |
(*_node_data)[vi].heap_index.erase(it); |
|
1040 |
if ((*_node_data)[vi].heap.empty()) { |
|
1041 |
_blossom_set->increase(v, std::numeric_limits<Value>::max()); |
|
1042 |
} else if ((*_blossom_set)[v] < (*_node_data)[vi].heap.prio()) { |
|
1043 |
_blossom_set->increase(v, (*_node_data)[vi].heap.prio()); |
|
1044 |
} |
|
1045 |
|
|
1046 |
if ((*_blossom_data)[vb].status == MATCHED) { |
|
1047 |
if (_blossom_set->classPrio(vb) == |
|
1048 |
std::numeric_limits<Value>::max()) { |
|
1049 |
_delta2->erase(vb); |
|
1050 |
} else if ((*_delta2)[vb] < _blossom_set->classPrio(vb) - |
|
1051 |
(*_blossom_data)[vb].offset) { |
|
1052 |
_delta2->increase(vb, _blossom_set->classPrio(vb) - |
|
1053 |
(*_blossom_data)[vb].offset); |
|
1054 |
} |
|
1055 |
} |
|
1186 | 1056 |
} |
1187 | 1057 |
} |
... | ... |
@@ -1190,9 +1060,38 @@ |
1190 | 1060 |
} |
1191 | 1061 |
|
1192 |
void |
|
1062 |
void oddToMatched(int blossom) { |
|
1063 |
(*_blossom_data)[blossom].offset -= _delta_sum; |
|
1064 |
|
|
1065 |
if (_blossom_set->classPrio(blossom) != |
|
1066 |
std::numeric_limits<Value>::max()) { |
|
1067 |
_delta2->push(blossom, _blossom_set->classPrio(blossom) - |
|
1068 |
(*_blossom_data)[blossom].offset); |
|
1069 |
} |
|
1070 |
|
|
1071 |
if (!_blossom_set->trivial(blossom)) { |
|
1072 |
_delta4->erase(blossom); |
|
1073 |
} |
|
1074 |
} |
|
1075 |
|
|
1076 |
void oddToEven(int blossom, int tree) { |
|
1077 |
if (!_blossom_set->trivial(blossom)) { |
|
1078 |
_delta4->erase(blossom); |
|
1079 |
(*_blossom_data)[blossom].pot -= |
|
1080 |
2 * (2 * _delta_sum - (*_blossom_data)[blossom].offset); |
|
1081 |
} |
|
1082 |
|
|
1193 | 1083 |
for (typename BlossomSet::ItemIt n(*_blossom_set, blossom); |
1194 | 1084 |
n != INVALID; ++n) { |
1195 | 1085 |
int ni = (*_node_index)[n]; |
1196 | 1086 |
|
1087 |
_blossom_set->increase(n, std::numeric_limits<Value>::max()); |
|
1088 |
|
|
1089 |
(*_node_data)[ni].heap.clear(); |
|
1090 |
(*_node_data)[ni].heap_index.clear(); |
|
1091 |
(*_node_data)[ni].pot += |
|
1092 |
2 * _delta_sum - (*_blossom_data)[blossom].offset; |
|
1093 |
|
|
1094 |
_delta1->push(n, (*_node_data)[ni].pot); |
|
1095 |
|
|
1197 | 1096 |
for (InArcIt e(_graph, n); e != INVALID; ++e) { |
1198 | 1097 |
Node v = _graph.source(e); |
... | ... |
@@ -1203,52 +1102,42 @@ |
1203 | 1102 |
dualScale * _weight[e]; |
1204 | 1103 |
|
1205 |
if (vb == blossom) { |
|
1206 |
if (_delta3->state(e) == _delta3->IN_HEAP) { |
|
1207 |
|
|
1104 |
if ((*_blossom_data)[vb].status == EVEN) { |
|
1105 |
if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) { |
|
1106 |
_delta3->push(e, rw / 2); |
|
1208 | 1107 |
} |
1209 |
} else if ((*_blossom_data)[vb].status == EVEN) { |
|
1210 |
|
|
1211 |
if (_delta3->state(e) == _delta3->IN_HEAP) { |
|
1212 |
_delta3->erase(e); |
|
1213 |
} |
|
1214 |
|
|
1215 |
int vt = _tree_set->find(vb); |
|
1216 |
|
|
1217 |
|
|
1108 |
} else { |
|
1218 | 1109 |
|
1219 | 1110 |
typename std::map<int, Arc>::iterator it = |
1220 |
(*_node_data)[ni].heap_index.find(vt); |
|
1221 |
|
|
1222 |
if (it != (*_node_data)[ni].heap_index.end()) { |
|
1223 |
if ((*_node_data)[ni].heap[it->second] > rw) { |
|
1224 |
(*_node_data)[ni].heap.replace(it->second, r); |
|
1225 |
(*_node_data)[ni].heap.decrease(r, rw); |
|
1226 |
|
|
1111 |
(*_node_data)[vi].heap_index.find(tree); |
|
1112 |
|
|
1113 |
if (it != (*_node_data)[vi].heap_index.end()) { |
|
1114 |
if ((*_node_data)[vi].heap[it->second] > rw) { |
|
1115 |
(*_node_data)[vi].heap.replace(it->second, e); |
|
1116 |
(*_node_data)[vi].heap.decrease(e, rw); |
|
1117 |
it->second = e; |
|
1227 | 1118 |
} |
1228 | 1119 |
} else { |
1229 |
(*_node_data)[ni].heap.push(r, rw); |
|
1230 |
(*_node_data)[ni].heap_index.insert(std::make_pair(vt, r)); |
|
1120 |
(*_node_data)[vi].heap.push(e, rw); |
|
1121 |
(*_node_data)[vi].heap_index.insert(std::make_pair(tree, e)); |
|
1231 | 1122 |
} |
1232 | 1123 |
|
1233 |
if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) { |
|
1234 |
_blossom_set->decrease(n, (*_node_data)[ni].heap.prio()); |
|
1235 |
|
|
1236 |
if (_delta2->state(blossom) != _delta2->IN_HEAP) { |
|
1237 |
_delta2->push(blossom, _blossom_set->classPrio(blossom) - |
|
1238 |
(*_blossom_data)[blossom].offset); |
|
1239 |
} else if ((*_delta2)[blossom] > _blossom_set->classPrio(blossom)- |
|
1240 |
(*_blossom_data)[blossom].offset){ |
|
1241 |
_delta2->decrease(blossom, _blossom_set->classPrio(blossom) - |
|
1242 |
(*_blossom_data)[blossom].offset); |
|
1124 |
if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) { |
|
1125 |
_blossom_set->decrease(v, (*_node_data)[vi].heap.prio()); |
|
1126 |
|
|
1127 |
if ((*_blossom_data)[vb].status == MATCHED) { |
|
1128 |
if (_delta2->state(vb) != _delta2->IN_HEAP) { |
|
1129 |
_delta2->push(vb, _blossom_set->classPrio(vb) - |
|
1130 |
(*_blossom_data)[vb].offset); |
|
1131 |
} else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) - |
|
1132 |
(*_blossom_data)[vb].offset) { |
|
1133 |
_delta2->decrease(vb, _blossom_set->classPrio(vb) - |
|
1134 |
(*_blossom_data)[vb].offset); |
|
1135 |
} |
|
1243 | 1136 |
} |
1244 | 1137 |
} |
1245 |
|
|
1246 |
} else if ((*_blossom_data)[vb].status == UNMATCHED) { |
|
1247 |
if (_delta3->state(e) == _delta3->IN_HEAP) { |
|
1248 |
_delta3->erase(e); |
|
1249 |
} |
|
1250 | 1138 |
} |
1251 | 1139 |
} |
1252 | 1140 |
} |
1141 |
(*_blossom_data)[blossom].offset = 0; |
|
1253 | 1142 |
} |
1254 | 1143 |
|
... | ... |
@@ -1295,10 +1184,8 @@ |
1295 | 1184 |
destroyTree(tree); |
1296 | 1185 |
|
1297 |
(*_blossom_data)[blossom].status = UNMATCHED; |
|
1298 | 1186 |
(*_blossom_data)[blossom].base = node; |
1299 |
|
|
1187 |
(*_blossom_data)[blossom].next = INVALID; |
|
1300 | 1188 |
} |
1301 | 1189 |
|
1302 |
|
|
1303 | 1190 |
void augmentOnEdge(const Edge& edge) { |
1304 | 1191 |
|
... | ... |
@@ -1306,21 +1193,11 @@ |
1306 | 1193 |
int right = _blossom_set->find(_graph.v(edge)); |
1307 | 1194 |
|
1308 |
if ((*_blossom_data)[left].status == EVEN) { |
|
1309 |
int left_tree = _tree_set->find(left); |
|
1310 |
alternatePath(left, left_tree); |
|
1311 |
destroyTree(left_tree); |
|
1312 |
} else { |
|
1313 |
(*_blossom_data)[left].status = MATCHED; |
|
1314 |
unmatchedToMatched(left); |
|
1315 |
} |
|
1316 |
|
|
1317 |
if ((*_blossom_data)[right].status == EVEN) { |
|
1318 |
int right_tree = _tree_set->find(right); |
|
1319 |
alternatePath(right, right_tree); |
|
1320 |
destroyTree(right_tree); |
|
1321 |
} else { |
|
1322 |
(*_blossom_data)[right].status = MATCHED; |
|
1323 |
unmatchedToMatched(right); |
|
1324 |
|
|
1195 |
int left_tree = _tree_set->find(left); |
|
1196 |
alternatePath(left, left_tree); |
|
1197 |
destroyTree(left_tree); |
|
1198 |
|
|
1199 |
int right_tree = _tree_set->find(right); |
|
1200 |
alternatePath(right, right_tree); |
|
1201 |
destroyTree(right_tree); |
|
1325 | 1202 |
|
1326 | 1203 |
(*_blossom_data)[left].next = _graph.direct(edge, true); |
... | ... |
@@ -1328,4 +1205,19 @@ |
1328 | 1205 |
} |
1329 | 1206 |
|
1207 |
void augmentOnArc(const Arc& arc) { |
|
1208 |
|
|
1209 |
int left = _blossom_set->find(_graph.source(arc)); |
|
1210 |
int right = _blossom_set->find(_graph.target(arc)); |
|
1211 |
|
|
1212 |
(*_blossom_data)[left].status = MATCHED; |
|
1213 |
|
|
1214 |
int right_tree = _tree_set->find(right); |
|
1215 |
alternatePath(right, right_tree); |
|
1216 |
destroyTree(right_tree); |
|
1217 |
|
|
1218 |
(*_blossom_data)[left].next = arc; |
|
1219 |
(*_blossom_data)[right].next = _graph.oppositeArc(arc); |
|
1220 |
} |
|
1221 |
|
|
1330 | 1222 |
void extendOnArc(const Arc& arc) { |
1331 | 1223 |
int base = _blossom_set->find(_graph.target(arc)); |
... | ... |
@@ -1530,5 +1422,5 @@ |
1530 | 1422 |
(*_blossom_data)[sb].pred = pred; |
1531 | 1423 |
(*_blossom_data)[sb].next = |
1532 |
|
|
1424 |
_graph.oppositeArc((*_blossom_data)[tb].next); |
|
1533 | 1425 |
|
1534 | 1426 |
pred = (*_blossom_data)[ub].next; |
... | ... |
@@ -1630,5 +1522,5 @@ |
1630 | 1522 |
|
1631 | 1523 |
for (int i = 0; i < int(blossoms.size()); ++i) { |
1632 |
if ((*_blossom_data)[blossoms[i]]. |
|
1524 |
if ((*_blossom_data)[blossoms[i]].next != INVALID) { |
|
1633 | 1525 |
|
1634 | 1526 |
Value offset = (*_blossom_data)[blossoms[i]].offset; |
... | ... |
@@ -1758,10 +1650,9 @@ |
1758 | 1650 |
_delta4->prio() : std::numeric_limits<Value>::max(); |
1759 | 1651 |
|
1760 |
_delta_sum = |
|
1652 |
_delta_sum = d3; OpType ot = D3; |
|
1653 |
if (d1 < _delta_sum) { _delta_sum = d1; ot = D1; } |
|
1761 | 1654 |
if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; } |
1762 |
if (d3 < _delta_sum) { _delta_sum = d3; ot = D3; } |
|
1763 | 1655 |
if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; } |
1764 | 1656 |
|
1765 |
|
|
1766 | 1657 |
switch (ot) { |
1767 | 1658 |
case D1: |
... | ... |
@@ -1776,6 +1667,11 @@ |
1776 | 1667 |
int blossom = _delta2->top(); |
1777 | 1668 |
Node n = _blossom_set->classTop(blossom); |
1778 |
Arc e = (*_node_data)[(*_node_index)[n]].heap.top(); |
|
1779 |
extendOnArc(e); |
|
1669 |
Arc a = (*_node_data)[(*_node_index)[n]].heap.top(); |
|
1670 |
if ((*_blossom_data)[blossom].next == INVALID) { |
|
1671 |
augmentOnArc(a); |
|
1672 |
--unmatched; |
|
1673 |
} else { |
|
1674 |
extendOnArc(a); |
|
1675 |
} |
|
1780 | 1676 |
} |
1781 | 1677 |
break; |
... | ... |
@@ -1790,18 +1686,6 @@ |
1790 | 1686 |
_delta3->pop(); |
1791 | 1687 |
} else { |
1792 |
int left_tree; |
|
1793 |
if ((*_blossom_data)[left_blossom].status == EVEN) { |
|
1794 |
left_tree = _tree_set->find(left_blossom); |
|
1795 |
} else { |
|
1796 |
left_tree = -1; |
|
1797 |
++unmatched; |
|
1798 |
} |
|
1799 |
int right_tree; |
|
1800 |
if ((*_blossom_data)[right_blossom].status == EVEN) { |
|
1801 |
right_tree = _tree_set->find(right_blossom); |
|
1802 |
} else { |
|
1803 |
right_tree = -1; |
|
1804 |
++unmatched; |
|
1805 |
} |
|
1688 |
int left_tree = _tree_set->find(left_blossom); |
|
1689 |
int right_tree = _tree_set->find(right_blossom); |
|
1806 | 1690 |
|
1807 | 1691 |
if (left_tree == right_tree) { |
... | ... |
@@ -1838,5 +1722,5 @@ |
1838 | 1722 |
|
1839 | 1723 |
/// \name Primal Solution |
1840 |
/// Functions to get the primal solution, i.e. the maximum weighted |
|
1724 |
/// Functions to get the primal solution, i.e. the maximum weighted |
|
1841 | 1725 |
/// matching.\n |
1842 | 1726 |
/// Either \ref run() or \ref start() function should be called before |
... | ... |
@@ -1857,5 +1741,5 @@ |
1857 | 1741 |
} |
1858 | 1742 |
} |
1859 |
return sum / |
|
1743 |
return sum / 2; |
|
1860 | 1744 |
} |
1861 | 1745 |
|
... | ... |
@@ -1877,5 +1761,5 @@ |
1877 | 1761 |
/// \brief Return \c true if the given edge is in the matching. |
1878 | 1762 |
/// |
1879 |
/// This function returns \c true if the given edge is in the found |
|
1763 |
/// This function returns \c true if the given edge is in the found |
|
1880 | 1764 |
/// matching. |
1881 | 1765 |
/// |
... | ... |
@@ -1888,5 +1772,5 @@ |
1888 | 1772 |
/// |
1889 | 1773 |
/// This function returns the matching arc (or edge) incident to the |
1890 |
/// given node in the found matching or \c INVALID if the node is |
|
1774 |
/// given node in the found matching or \c INVALID if the node is |
|
1891 | 1775 |
/// not covered by the matching. |
1892 | 1776 |
/// |
... | ... |
@@ -1906,5 +1790,5 @@ |
1906 | 1790 |
/// \brief Return the mate of the given node. |
1907 | 1791 |
/// |
1908 |
/// This function returns the mate of the given node in the found |
|
1792 |
/// This function returns the mate of the given node in the found |
|
1909 | 1793 |
/// matching or \c INVALID if the node is not covered by the matching. |
1910 | 1794 |
/// |
... | ... |
@@ -1926,6 +1810,6 @@ |
1926 | 1810 |
/// \brief Return the value of the dual solution. |
1927 | 1811 |
/// |
1928 |
/// This function returns the value of the dual solution. |
|
1929 |
/// It should be equal to the primal value scaled by \ref dualScale |
|
1812 |
/// This function returns the value of the dual solution. |
|
1813 |
/// It should be equal to the primal value scaled by \ref dualScale |
|
1930 | 1814 |
/// "dual scale". |
1931 | 1815 |
/// |
... | ... |
@@ -1982,7 +1866,7 @@ |
1982 | 1866 |
/// \brief Iterator for obtaining the nodes of a blossom. |
1983 | 1867 |
/// |
1984 |
/// This class provides an iterator for obtaining the nodes of the |
|
1868 |
/// This class provides an iterator for obtaining the nodes of the |
|
1985 | 1869 |
/// given blossom. It lists a subset of the nodes. |
1986 |
/// Before using this iterator, you must allocate a |
|
1870 |
/// Before using this iterator, you must allocate a |
|
1987 | 1871 |
/// MaxWeightedMatching class and execute it. |
1988 | 1872 |
class BlossomIt { |
... | ... |
@@ -1993,6 +1877,6 @@ |
1993 | 1877 |
/// Constructor to get the nodes of the given variable. |
1994 | 1878 |
/// |
1995 |
/// \pre Either \ref MaxWeightedMatching::run() "algorithm.run()" or |
|
1996 |
/// \ref MaxWeightedMatching::start() "algorithm.start()" must be |
|
1879 |
/// \pre Either \ref MaxWeightedMatching::run() "algorithm.run()" or |
|
1880 |
/// \ref MaxWeightedMatching::start() "algorithm.start()" must be |
|
1997 | 1881 |
/// called before initializing this iterator. |
1998 | 1882 |
BlossomIt(const MaxWeightedMatching& algorithm, int variable) |
... | ... |
@@ -2047,6 +1931,6 @@ |
2047 | 1931 |
/// \f$O(nm\log n)\f$ time complexity. |
2048 | 1932 |
/// |
2049 |
/// The maximum weighted perfect matching problem is to find a subset of |
|
2050 |
/// the edges in an undirected graph with maximum overall weight for which |
|
1933 |
/// The maximum weighted perfect matching problem is to find a subset of |
|
1934 |
/// the edges in an undirected graph with maximum overall weight for which |
|
2051 | 1935 |
/// each node has exactly one incident edge. |
2052 | 1936 |
/// It can be formulated with the following linear program. |
... | ... |
@@ -2071,14 +1955,14 @@ |
2071 | 1955 |
\frac{\vert B \vert - 1}{2}z_B\f] */ |
2072 | 1956 |
/// |
2073 |
/// The algorithm can be executed with the run() function. |
|
1957 |
/// The algorithm can be executed with the run() function. |
|
2074 | 1958 |
/// After it the matching (the primal solution) and the dual solution |
2075 |
/// can be obtained using the query functions and the |
|
2076 |
/// \ref MaxWeightedPerfectMatching::BlossomIt "BlossomIt" nested class, |
|
2077 |
/// |
|
1959 |
/// can be obtained using the query functions and the |
|
1960 |
/// \ref MaxWeightedPerfectMatching::BlossomIt "BlossomIt" nested class, |
|
1961 |
/// which is able to iterate on the nodes of a blossom. |
|
2078 | 1962 |
/// If the value type is integer, then the dual solution is multiplied |
2079 | 1963 |
/// by \ref MaxWeightedMatching::dualScale "4". |
2080 | 1964 |
/// |
2081 | 1965 |
/// \tparam GR The undirected graph type the algorithm runs on. |
2082 |
/// \tparam WM The type edge weight map. The default type is |
|
1966 |
/// \tparam WM The type edge weight map. The default type is |
|
2083 | 1967 |
/// \ref concepts::Graph::EdgeMap "GR::EdgeMap<int>". |
2084 | 1968 |
#ifdef DOXYGEN |
... | ... |
@@ -2234,7 +2118,4 @@ |
2234 | 2118 |
|
2235 | 2119 |
void destroyStructures() { |
2236 |
_node_num = countNodes(_graph); |
|
2237 |
_blossom_num = _node_num * 3 / 2; |
|
2238 |
|
|
2239 | 2120 |
if (_matching) { |
2240 | 2121 |
delete _matching; |
... | ... |
@@ -2992,6 +2873,6 @@ |
2992 | 2873 |
_delta4->prio() : std::numeric_limits<Value>::max(); |
2993 | 2874 |
|
2994 |
_delta_sum = d2; OpType ot = D2; |
|
2995 |
if (d3 < _delta_sum) { _delta_sum = d3; ot = D3; } |
|
2875 |
_delta_sum = d3; OpType ot = D3; |
|
2876 |
if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; } |
|
2996 | 2877 |
if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; } |
2997 | 2878 |
|
... | ... |
@@ -3056,5 +2937,5 @@ |
3056 | 2937 |
|
3057 | 2938 |
/// \name Primal Solution |
3058 |
/// Functions to get the primal solution, i.e. the maximum weighted |
|
2939 |
/// Functions to get the primal solution, i.e. the maximum weighted |
|
3059 | 2940 |
/// perfect matching.\n |
3060 | 2941 |
/// Either \ref run() or \ref start() function should be called before |
... | ... |
@@ -3075,10 +2956,10 @@ |
3075 | 2956 |
} |
3076 | 2957 |
} |
3077 |
return sum / |
|
2958 |
return sum / 2; |
|
3078 | 2959 |
} |
3079 | 2960 |
|
3080 | 2961 |
/// \brief Return \c true if the given edge is in the matching. |
3081 | 2962 |
/// |
3082 |
/// This function returns \c true if the given edge is in the found |
|
2963 |
/// This function returns \c true if the given edge is in the found |
|
3083 | 2964 |
/// matching. |
3084 | 2965 |
/// |
... | ... |
@@ -3091,5 +2972,5 @@ |
3091 | 2972 |
/// |
3092 | 2973 |
/// This function returns the matching arc (or edge) incident to the |
3093 |
/// given node in the found matching or \c INVALID if the node is |
|
2974 |
/// given node in the found matching or \c INVALID if the node is |
|
3094 | 2975 |
/// not covered by the matching. |
3095 | 2976 |
/// |
... | ... |
@@ -3109,5 +2990,5 @@ |
3109 | 2990 |
/// \brief Return the mate of the given node. |
3110 | 2991 |
/// |
3111 |
/// This function returns the mate of the given node in the found |
|
2992 |
/// This function returns the mate of the given node in the found |
|
3112 | 2993 |
/// matching or \c INVALID if the node is not covered by the matching. |
3113 | 2994 |
/// |
... | ... |
@@ -3128,6 +3009,6 @@ |
3128 | 3009 |
/// \brief Return the value of the dual solution. |
3129 | 3010 |
/// |
3130 |
/// This function returns the value of the dual solution. |
|
3131 |
/// It should be equal to the primal value scaled by \ref dualScale |
|
3011 |
/// This function returns the value of the dual solution. |
|
3012 |
/// It should be equal to the primal value scaled by \ref dualScale |
|
3132 | 3013 |
/// "dual scale". |
3133 | 3014 |
/// |
... | ... |
@@ -3184,7 +3065,7 @@ |
3184 | 3065 |
/// \brief Iterator for obtaining the nodes of a blossom. |
3185 | 3066 |
/// |
3186 |
/// This class provides an iterator for obtaining the nodes of the |
|
3067 |
/// This class provides an iterator for obtaining the nodes of the |
|
3187 | 3068 |
/// given blossom. It lists a subset of the nodes. |
3188 |
/// Before using this iterator, you must allocate a |
|
3069 |
/// Before using this iterator, you must allocate a |
|
3189 | 3070 |
/// MaxWeightedPerfectMatching class and execute it. |
3190 | 3071 |
class BlossomIt { |
... | ... |
@@ -3195,6 +3076,6 @@ |
3195 | 3076 |
/// Constructor to get the nodes of the given variable. |
3196 | 3077 |
/// |
3197 |
/// \pre Either \ref MaxWeightedPerfectMatching::run() "algorithm.run()" |
|
3198 |
/// or \ref MaxWeightedPerfectMatching::start() "algorithm.start()" |
|
3078 |
/// \pre Either \ref MaxWeightedPerfectMatching::run() "algorithm.run()" |
|
3079 |
/// or \ref MaxWeightedPerfectMatching::start() "algorithm.start()" |
|
3199 | 3080 |
/// must be called before initializing this iterator. |
3200 | 3081 |
BlossomIt(const MaxWeightedPerfectMatching& algorithm, int variable) |
... | ... |
@@ -3242,3 +3123,3 @@ |
3242 | 3123 |
} //END OF NAMESPACE LEMON |
3243 | 3124 |
|
3244 |
#endif // |
|
3125 |
#endif //LEMON_MATCHING_H |
0 comments (0 inline)