39 /// |
39 /// |
40 /// \brief Maximum cardinality matching in general graphs |
40 /// \brief Maximum cardinality matching in general graphs |
41 /// |
41 /// |
42 /// This class implements Edmonds' alternating forest matching algorithm |
42 /// This class implements Edmonds' alternating forest matching algorithm |
43 /// for finding a maximum cardinality matching in a general undirected graph. |
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 /// (the default is the empty one). |
45 /// (the default is the empty one). |
46 /// |
46 /// |
47 /// The dual solution of the problem is a map of the nodes to |
47 /// The dual solution of the problem is a map of the nodes to |
48 /// \ref MaxMatching::Status "Status", having values \c EVEN (or \c D), |
48 /// \ref MaxMatching::Status "Status", having values \c EVEN (or \c D), |
49 /// \c ODD (or \c A) and \c MATCHED (or \c C) defining the Gallai-Edmonds |
49 /// \c ODD (or \c A) and \c MATCHED (or \c C) defining the Gallai-Edmonds |
67 typedef typename Graph::template NodeMap<typename Graph::Arc> |
67 typedef typename Graph::template NodeMap<typename Graph::Arc> |
68 MatchingMap; |
68 MatchingMap; |
69 |
69 |
70 ///\brief Status constants for Gallai-Edmonds decomposition. |
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 ///decomposition of a graph. The nodes with status \c EVEN (or \c D) |
73 ///decomposition of a graph. The nodes with status \c EVEN (or \c D) |
74 ///induce a subgraph with factor-critical components, the nodes with |
74 ///induce a subgraph with factor-critical components, the nodes with |
75 ///status \c ODD (or \c A) form the canonical barrier, and the nodes |
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 ///perfect matching. |
77 ///perfect matching. |
78 enum Status { |
78 enum Status { |
79 EVEN = 1, ///< = 1. (\c D is an alias for \c EVEN.) |
79 EVEN = 1, ///< = 1. (\c D is an alias for \c EVEN.) |
80 D = 1, |
80 D = 1, |
81 MATCHED = 0, ///< = 0. (\c C is an alias for \c MATCHED.) |
81 MATCHED = 0, ///< = 0. (\c C is an alias for \c MATCHED.) |
568 return size / 2; |
568 return size / 2; |
569 } |
569 } |
570 |
570 |
571 /// \brief Return \c true if the given edge is in the matching. |
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 /// matching. |
574 /// matching. |
575 bool matching(const Edge& edge) const { |
575 bool matching(const Edge& edge) const { |
576 return edge == (*_matching)[_graph.u(edge)]; |
576 return edge == (*_matching)[_graph.u(edge)]; |
577 } |
577 } |
578 |
578 |
579 /// \brief Return the matching arc (or edge) incident to the given node. |
579 /// \brief Return the matching arc (or edge) incident to the given node. |
580 /// |
580 /// |
581 /// This function returns the matching arc (or edge) incident to the |
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 /// not covered by the matching. |
583 /// not covered by the matching. |
584 Arc matching(const Node& n) const { |
584 Arc matching(const Node& n) const { |
585 return (*_matching)[n]; |
585 return (*_matching)[n]; |
586 } |
586 } |
587 |
587 |
593 return *_matching; |
593 return *_matching; |
594 } |
594 } |
595 |
595 |
596 /// \brief Return the mate of the given node. |
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 /// matching or \c INVALID if the node is not covered by the matching. |
599 /// matching or \c INVALID if the node is not covered by the matching. |
600 Node mate(const Node& n) const { |
600 Node mate(const Node& n) const { |
601 return (*_matching)[n] != INVALID ? |
601 return (*_matching)[n] != INVALID ? |
602 _graph.target((*_matching)[n]) : INVALID; |
602 _graph.target((*_matching)[n]) : INVALID; |
603 } |
603 } |
604 |
604 |
605 /// @} |
605 /// @} |
606 |
606 |
607 /// \name Dual Solution |
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 /// decomposition. |
609 /// decomposition. |
610 |
610 |
611 /// @{ |
611 /// @{ |
612 |
612 |
613 /// \brief Return the status of the given node in the Edmonds-Gallai |
613 /// \brief Return the status of the given node in the Edmonds-Gallai |
646 /// This class provides an efficient implementation of Edmond's |
646 /// This class provides an efficient implementation of Edmond's |
647 /// maximum weighted matching algorithm. The implementation is based |
647 /// maximum weighted matching algorithm. The implementation is based |
648 /// on extensive use of priority queues and provides |
648 /// on extensive use of priority queues and provides |
649 /// \f$O(nm\log n)\f$ time complexity. |
649 /// \f$O(nm\log n)\f$ time complexity. |
650 /// |
650 /// |
651 /// The maximum weighted matching problem is to find a subset of the |
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 |
652 /// edges in an undirected graph with maximum overall weight for which |
653 /// each node has at most one incident edge. |
653 /// each node has at most one incident edge. |
654 /// It can be formulated with the following linear program. |
654 /// It can be formulated with the following linear program. |
655 /// \f[ \sum_{e \in \delta(u)}x_e \le 1 \quad \forall u\in V\f] |
655 /// \f[ \sum_{e \in \delta(u)}x_e \le 1 \quad \forall u\in V\f] |
656 /** \f[ \sum_{e \in \gamma(B)}x_e \le \frac{\vert B \vert - 1}{2} |
656 /** \f[ \sum_{e \in \gamma(B)}x_e \le \frac{\vert B \vert - 1}{2} |
657 \quad \forall B\in\mathcal{O}\f] */ |
657 \quad \forall B\in\mathcal{O}\f] */ |
671 /// \f[y_u \ge 0 \quad \forall u \in V\f] |
671 /// \f[y_u \ge 0 \quad \forall u \in V\f] |
672 /// \f[z_B \ge 0 \quad \forall B \in \mathcal{O}\f] |
672 /// \f[z_B \ge 0 \quad \forall B \in \mathcal{O}\f] |
673 /** \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}} |
673 /** \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}} |
674 \frac{\vert B \vert - 1}{2}z_B\f] */ |
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 /// After it the matching (the primal solution) and the dual solution |
677 /// After it the matching (the primal solution) and the dual solution |
678 /// can be obtained using the query functions and the |
678 /// can be obtained using the query functions and the |
679 /// \ref MaxWeightedMatching::BlossomIt "BlossomIt" nested class, |
679 /// \ref MaxWeightedMatching::BlossomIt "BlossomIt" nested class, |
680 /// which is able to iterate on the nodes of a blossom. |
680 /// which is able to iterate on the nodes of a blossom. |
681 /// If the value type is integer, then the dual solution is multiplied |
681 /// If the value type is integer, then the dual solution is multiplied |
682 /// by \ref MaxWeightedMatching::dualScale "4". |
682 /// by \ref MaxWeightedMatching::dualScale "4". |
683 /// |
683 /// |
684 /// \tparam GR The undirected graph type the algorithm runs on. |
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 /// \ref concepts::Graph::EdgeMap "GR::EdgeMap<int>". |
686 /// \ref concepts::Graph::EdgeMap "GR::EdgeMap<int>". |
687 #ifdef DOXYGEN |
687 #ifdef DOXYGEN |
688 template <typename GR, typename WM> |
688 template <typename GR, typename WM> |
689 #else |
689 #else |
690 template <typename GR, |
690 template <typename GR, |
919 dualScale * _weight[e]; |
916 dualScale * _weight[e]; |
920 |
917 |
921 if ((*_blossom_data)[vb].status == EVEN) { |
918 if ((*_blossom_data)[vb].status == EVEN) { |
922 if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) { |
919 if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) { |
923 _delta3->push(e, rw / 2); |
920 _delta3->push(e, rw / 2); |
924 } |
|
925 } else if ((*_blossom_data)[vb].status == UNMATCHED) { |
|
926 if (_delta3->state(e) != _delta3->IN_HEAP) { |
|
927 _delta3->push(e, rw); |
|
928 } |
921 } |
929 } else { |
922 } else { |
930 typename std::map<int, Arc>::iterator it = |
923 typename std::map<int, Arc>::iterator it = |
931 (*_node_data)[vi].heap_index.find(tree); |
924 (*_node_data)[vi].heap_index.find(tree); |
932 |
925 |
947 if ((*_blossom_data)[vb].status == MATCHED) { |
940 if ((*_blossom_data)[vb].status == MATCHED) { |
948 if (_delta2->state(vb) != _delta2->IN_HEAP) { |
941 if (_delta2->state(vb) != _delta2->IN_HEAP) { |
949 _delta2->push(vb, _blossom_set->classPrio(vb) - |
942 _delta2->push(vb, _blossom_set->classPrio(vb) - |
950 (*_blossom_data)[vb].offset); |
943 (*_blossom_data)[vb].offset); |
951 } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) - |
944 } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) - |
952 (*_blossom_data)[vb].offset){ |
945 (*_blossom_data)[vb].offset) { |
953 _delta2->decrease(vb, _blossom_set->classPrio(vb) - |
946 _delta2->decrease(vb, _blossom_set->classPrio(vb) - |
954 (*_blossom_data)[vb].offset); |
947 (*_blossom_data)[vb].offset); |
955 } |
948 } |
956 } |
949 } |
957 } |
950 } |
966 _delta2->erase(blossom); |
959 _delta2->erase(blossom); |
967 } |
960 } |
968 (*_blossom_data)[blossom].offset += _delta_sum; |
961 (*_blossom_data)[blossom].offset += _delta_sum; |
969 if (!_blossom_set->trivial(blossom)) { |
962 if (!_blossom_set->trivial(blossom)) { |
970 _delta4->push(blossom, (*_blossom_data)[blossom].pot / 2 + |
963 _delta4->push(blossom, (*_blossom_data)[blossom].pot / 2 + |
971 (*_blossom_data)[blossom].offset); |
964 (*_blossom_data)[blossom].offset); |
972 } |
965 } |
973 } |
966 } |
974 |
967 |
975 void evenToMatched(int blossom, int tree) { |
968 void evenToMatched(int blossom, int tree) { |
976 if (!_blossom_set->trivial(blossom)) { |
969 if (!_blossom_set->trivial(blossom)) { |
1075 (*_blossom_data)[blossom].offset -= _delta_sum; |
1063 (*_blossom_data)[blossom].offset -= _delta_sum; |
1076 |
1064 |
1077 if (_blossom_set->classPrio(blossom) != |
1065 if (_blossom_set->classPrio(blossom) != |
1078 std::numeric_limits<Value>::max()) { |
1066 std::numeric_limits<Value>::max()) { |
1079 _delta2->push(blossom, _blossom_set->classPrio(blossom) - |
1067 _delta2->push(blossom, _blossom_set->classPrio(blossom) - |
1080 (*_blossom_data)[blossom].offset); |
1068 (*_blossom_data)[blossom].offset); |
1081 } |
1069 } |
1082 |
1070 |
1083 if (!_blossom_set->trivial(blossom)) { |
1071 if (!_blossom_set->trivial(blossom)) { |
1084 _delta4->erase(blossom); |
1072 _delta4->erase(blossom); |
1085 } |
1073 } |
1114 dualScale * _weight[e]; |
1102 dualScale * _weight[e]; |
1115 |
1103 |
1116 if ((*_blossom_data)[vb].status == EVEN) { |
1104 if ((*_blossom_data)[vb].status == EVEN) { |
1117 if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) { |
1105 if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) { |
1118 _delta3->push(e, rw / 2); |
1106 _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 } |
1107 } |
1124 } else { |
1108 } else { |
1125 |
1109 |
1126 typename std::map<int, Arc>::iterator it = |
1110 typename std::map<int, Arc>::iterator it = |
1127 (*_node_data)[vi].heap_index.find(tree); |
1111 (*_node_data)[vi].heap_index.find(tree); |
1155 } |
1139 } |
1156 } |
1140 } |
1157 (*_blossom_data)[blossom].offset = 0; |
1141 (*_blossom_data)[blossom].offset = 0; |
1158 } |
1142 } |
1159 |
1143 |
1160 |
|
1161 void matchedToUnmatched(int blossom) { |
|
1162 if (_delta2->state(blossom) == _delta2->IN_HEAP) { |
|
1163 _delta2->erase(blossom); |
|
1164 } |
|
1165 |
|
1166 for (typename BlossomSet::ItemIt n(*_blossom_set, blossom); |
|
1167 n != INVALID; ++n) { |
|
1168 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); |
|
1177 int vb = _blossom_set->find(v); |
|
1178 int vi = (*_node_index)[v]; |
|
1179 |
|
1180 Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot - |
|
1181 dualScale * _weight[e]; |
|
1182 |
|
1183 if ((*_blossom_data)[vb].status == EVEN) { |
|
1184 if (_delta3->state(e) != _delta3->IN_HEAP) { |
|
1185 _delta3->push(e, rw); |
|
1186 } |
|
1187 } |
|
1188 } |
|
1189 } |
|
1190 } |
|
1191 |
|
1192 void unmatchedToMatched(int blossom) { |
|
1193 for (typename BlossomSet::ItemIt n(*_blossom_set, blossom); |
|
1194 n != INVALID; ++n) { |
|
1195 int ni = (*_node_index)[n]; |
|
1196 |
|
1197 for (InArcIt e(_graph, n); e != INVALID; ++e) { |
|
1198 Node v = _graph.source(e); |
|
1199 int vb = _blossom_set->find(v); |
|
1200 int vi = (*_node_index)[v]; |
|
1201 |
|
1202 Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot - |
|
1203 dualScale * _weight[e]; |
|
1204 |
|
1205 if (vb == blossom) { |
|
1206 if (_delta3->state(e) == _delta3->IN_HEAP) { |
|
1207 _delta3->erase(e); |
|
1208 } |
|
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 Arc r = _graph.oppositeArc(e); |
|
1218 |
|
1219 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 it->second = r; |
|
1227 } |
|
1228 } else { |
|
1229 (*_node_data)[ni].heap.push(r, rw); |
|
1230 (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r)); |
|
1231 } |
|
1232 |
|
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); |
|
1243 } |
|
1244 } |
|
1245 |
|
1246 } else if ((*_blossom_data)[vb].status == UNMATCHED) { |
|
1247 if (_delta3->state(e) == _delta3->IN_HEAP) { |
|
1248 _delta3->erase(e); |
|
1249 } |
|
1250 } |
|
1251 } |
|
1252 } |
|
1253 } |
|
1254 |
|
1255 void alternatePath(int even, int tree) { |
1144 void alternatePath(int even, int tree) { |
1256 int odd; |
1145 int odd; |
1257 |
1146 |
1258 evenToMatched(even, tree); |
1147 evenToMatched(even, tree); |
1259 (*_blossom_data)[even].status = MATCHED; |
1148 (*_blossom_data)[even].status = MATCHED; |
1292 int tree = _tree_set->find(blossom); |
1181 int tree = _tree_set->find(blossom); |
1293 |
1182 |
1294 alternatePath(blossom, tree); |
1183 alternatePath(blossom, tree); |
1295 destroyTree(tree); |
1184 destroyTree(tree); |
1296 |
1185 |
1297 (*_blossom_data)[blossom].status = UNMATCHED; |
|
1298 (*_blossom_data)[blossom].base = node; |
1186 (*_blossom_data)[blossom].base = node; |
1299 matchedToUnmatched(blossom); |
1187 (*_blossom_data)[blossom].next = INVALID; |
1300 } |
1188 } |
1301 |
|
1302 |
1189 |
1303 void augmentOnEdge(const Edge& edge) { |
1190 void augmentOnEdge(const Edge& edge) { |
1304 |
1191 |
1305 int left = _blossom_set->find(_graph.u(edge)); |
1192 int left = _blossom_set->find(_graph.u(edge)); |
1306 int right = _blossom_set->find(_graph.v(edge)); |
1193 int right = _blossom_set->find(_graph.v(edge)); |
1307 |
1194 |
1308 if ((*_blossom_data)[left].status == EVEN) { |
1195 int left_tree = _tree_set->find(left); |
1309 int left_tree = _tree_set->find(left); |
1196 alternatePath(left, left_tree); |
1310 alternatePath(left, left_tree); |
1197 destroyTree(left_tree); |
1311 destroyTree(left_tree); |
1198 |
1312 } else { |
1199 int right_tree = _tree_set->find(right); |
1313 (*_blossom_data)[left].status = MATCHED; |
1200 alternatePath(right, right_tree); |
1314 unmatchedToMatched(left); |
1201 destroyTree(right_tree); |
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 } |
|
1325 |
1202 |
1326 (*_blossom_data)[left].next = _graph.direct(edge, true); |
1203 (*_blossom_data)[left].next = _graph.direct(edge, true); |
1327 (*_blossom_data)[right].next = _graph.direct(edge, false); |
1204 (*_blossom_data)[right].next = _graph.direct(edge, false); |
|
1205 } |
|
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); |
1328 } |
1220 } |
1329 |
1221 |
1330 void extendOnArc(const Arc& arc) { |
1222 void extendOnArc(const Arc& arc) { |
1331 int base = _blossom_set->find(_graph.target(arc)); |
1223 int base = _blossom_set->find(_graph.target(arc)); |
1332 int tree = _tree_set->find(base); |
1224 int tree = _tree_set->find(base); |
1527 (*_blossom_data)[sb].status = ODD; |
1419 (*_blossom_data)[sb].status = ODD; |
1528 matchedToOdd(sb); |
1420 matchedToOdd(sb); |
1529 _tree_set->insert(sb, tree); |
1421 _tree_set->insert(sb, tree); |
1530 (*_blossom_data)[sb].pred = pred; |
1422 (*_blossom_data)[sb].pred = pred; |
1531 (*_blossom_data)[sb].next = |
1423 (*_blossom_data)[sb].next = |
1532 _graph.oppositeArc((*_blossom_data)[tb].next); |
1424 _graph.oppositeArc((*_blossom_data)[tb].next); |
1533 |
1425 |
1534 pred = (*_blossom_data)[ub].next; |
1426 pred = (*_blossom_data)[ub].next; |
1535 |
1427 |
1536 (*_blossom_data)[tb].status = EVEN; |
1428 (*_blossom_data)[tb].status = EVEN; |
1537 matchedToEven(tb, tree); |
1429 matchedToEven(tb, tree); |
1627 for (typename BlossomSet::ClassIt c(*_blossom_set); c != INVALID; ++c) { |
1519 for (typename BlossomSet::ClassIt c(*_blossom_set); c != INVALID; ++c) { |
1628 blossoms.push_back(c); |
1520 blossoms.push_back(c); |
1629 } |
1521 } |
1630 |
1522 |
1631 for (int i = 0; i < int(blossoms.size()); ++i) { |
1523 for (int i = 0; i < int(blossoms.size()); ++i) { |
1632 if ((*_blossom_data)[blossoms[i]].status == MATCHED) { |
1524 if ((*_blossom_data)[blossoms[i]].next != INVALID) { |
1633 |
1525 |
1634 Value offset = (*_blossom_data)[blossoms[i]].offset; |
1526 Value offset = (*_blossom_data)[blossoms[i]].offset; |
1635 (*_blossom_data)[blossoms[i]].pot += 2 * offset; |
1527 (*_blossom_data)[blossoms[i]].pot += 2 * offset; |
1636 for (typename BlossomSet::ItemIt n(*_blossom_set, blossoms[i]); |
1528 for (typename BlossomSet::ItemIt n(*_blossom_set, blossoms[i]); |
1637 n != INVALID; ++n) { |
1529 n != INVALID; ++n) { |
1755 _delta3->prio() : std::numeric_limits<Value>::max(); |
1647 _delta3->prio() : std::numeric_limits<Value>::max(); |
1756 |
1648 |
1757 Value d4 = !_delta4->empty() ? |
1649 Value d4 = !_delta4->empty() ? |
1758 _delta4->prio() : std::numeric_limits<Value>::max(); |
1650 _delta4->prio() : std::numeric_limits<Value>::max(); |
1759 |
1651 |
1760 _delta_sum = d1; OpType ot = D1; |
1652 _delta_sum = d3; OpType ot = D3; |
|
1653 if (d1 < _delta_sum) { _delta_sum = d1; ot = D1; } |
1761 if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; } |
1654 if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; } |
1762 if (d3 < _delta_sum) { _delta_sum = d3; ot = D3; } |
|
1763 if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; } |
1655 if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; } |
1764 |
|
1765 |
1656 |
1766 switch (ot) { |
1657 switch (ot) { |
1767 case D1: |
1658 case D1: |
1768 { |
1659 { |
1769 Node n = _delta1->top(); |
1660 Node n = _delta1->top(); |
1787 int right_blossom = _blossom_set->find(_graph.v(e)); |
1683 int right_blossom = _blossom_set->find(_graph.v(e)); |
1788 |
1684 |
1789 if (left_blossom == right_blossom) { |
1685 if (left_blossom == right_blossom) { |
1790 _delta3->pop(); |
1686 _delta3->pop(); |
1791 } else { |
1687 } else { |
1792 int left_tree; |
1688 int left_tree = _tree_set->find(left_blossom); |
1793 if ((*_blossom_data)[left_blossom].status == EVEN) { |
1689 int right_tree = _tree_set->find(right_blossom); |
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 } |
|
1806 |
1690 |
1807 if (left_tree == right_tree) { |
1691 if (left_tree == right_tree) { |
1808 shrinkOnEdge(e, left_tree); |
1692 shrinkOnEdge(e, left_tree); |
1809 } else { |
1693 } else { |
1810 augmentOnEdge(e); |
1694 augmentOnEdge(e); |
1874 return num /= 2; |
1758 return num /= 2; |
1875 } |
1759 } |
1876 |
1760 |
1877 /// \brief Return \c true if the given edge is in the matching. |
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 /// matching. |
1764 /// matching. |
1881 /// |
1765 /// |
1882 /// \pre Either run() or start() must be called before using this function. |
1766 /// \pre Either run() or start() must be called before using this function. |
1883 bool matching(const Edge& edge) const { |
1767 bool matching(const Edge& edge) const { |
1884 return edge == (*_matching)[_graph.u(edge)]; |
1768 return edge == (*_matching)[_graph.u(edge)]; |
1885 } |
1769 } |
1886 |
1770 |
1887 /// \brief Return the matching arc (or edge) incident to the given node. |
1771 /// \brief Return the matching arc (or edge) incident to the given node. |
1888 /// |
1772 /// |
1889 /// This function returns the matching arc (or edge) incident to the |
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 /// not covered by the matching. |
1775 /// not covered by the matching. |
1892 /// |
1776 /// |
1893 /// \pre Either run() or start() must be called before using this function. |
1777 /// \pre Either run() or start() must be called before using this function. |
1894 Arc matching(const Node& node) const { |
1778 Arc matching(const Node& node) const { |
1895 return (*_matching)[node]; |
1779 return (*_matching)[node]; |
1979 return _blossom_potential[k].value; |
1863 return _blossom_potential[k].value; |
1980 } |
1864 } |
1981 |
1865 |
1982 /// \brief Iterator for obtaining the nodes of a blossom. |
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 /// given blossom. It lists a subset of the nodes. |
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 /// MaxWeightedMatching class and execute it. |
1871 /// MaxWeightedMatching class and execute it. |
1988 class BlossomIt { |
1872 class BlossomIt { |
1989 public: |
1873 public: |
1990 |
1874 |
1991 /// \brief Constructor. |
1875 /// \brief Constructor. |
1992 /// |
1876 /// |
1993 /// Constructor to get the nodes of the given variable. |
1877 /// Constructor to get the nodes of the given variable. |
1994 /// |
1878 /// |
1995 /// \pre Either \ref MaxWeightedMatching::run() "algorithm.run()" or |
1879 /// \pre Either \ref MaxWeightedMatching::run() "algorithm.run()" or |
1996 /// \ref MaxWeightedMatching::start() "algorithm.start()" must be |
1880 /// \ref MaxWeightedMatching::start() "algorithm.start()" must be |
1997 /// called before initializing this iterator. |
1881 /// called before initializing this iterator. |
1998 BlossomIt(const MaxWeightedMatching& algorithm, int variable) |
1882 BlossomIt(const MaxWeightedMatching& algorithm, int variable) |
1999 : _algorithm(&algorithm) |
1883 : _algorithm(&algorithm) |
2000 { |
1884 { |
2001 _index = _algorithm->_blossom_potential[variable].begin; |
1885 _index = _algorithm->_blossom_potential[variable].begin; |
2044 /// This class provides an efficient implementation of Edmond's |
1928 /// This class provides an efficient implementation of Edmond's |
2045 /// maximum weighted perfect matching algorithm. The implementation |
1929 /// maximum weighted perfect matching algorithm. The implementation |
2046 /// is based on extensive use of priority queues and provides |
1930 /// is based on extensive use of priority queues and provides |
2047 /// \f$O(nm\log n)\f$ time complexity. |
1931 /// \f$O(nm\log n)\f$ time complexity. |
2048 /// |
1932 /// |
2049 /// The maximum weighted perfect matching problem is to find a subset of |
1933 /// 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 |
1934 /// the edges in an undirected graph with maximum overall weight for which |
2051 /// each node has exactly one incident edge. |
1935 /// each node has exactly one incident edge. |
2052 /// It can be formulated with the following linear program. |
1936 /// It can be formulated with the following linear program. |
2053 /// \f[ \sum_{e \in \delta(u)}x_e = 1 \quad \forall u\in V\f] |
1937 /// \f[ \sum_{e \in \delta(u)}x_e = 1 \quad \forall u\in V\f] |
2054 /** \f[ \sum_{e \in \gamma(B)}x_e \le \frac{\vert B \vert - 1}{2} |
1938 /** \f[ \sum_{e \in \gamma(B)}x_e \le \frac{\vert B \vert - 1}{2} |
2055 \quad \forall B\in\mathcal{O}\f] */ |
1939 \quad \forall B\in\mathcal{O}\f] */ |
2068 w_{uv} \quad \forall uv\in E\f] */ |
1952 w_{uv} \quad \forall uv\in E\f] */ |
2069 /// \f[z_B \ge 0 \quad \forall B \in \mathcal{O}\f] |
1953 /// \f[z_B \ge 0 \quad \forall B \in \mathcal{O}\f] |
2070 /** \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}} |
1954 /** \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}} |
2071 \frac{\vert B \vert - 1}{2}z_B\f] */ |
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 /// After it the matching (the primal solution) and the dual solution |
1958 /// After it the matching (the primal solution) and the dual solution |
2075 /// can be obtained using the query functions and the |
1959 /// can be obtained using the query functions and the |
2076 /// \ref MaxWeightedPerfectMatching::BlossomIt "BlossomIt" nested class, |
1960 /// \ref MaxWeightedPerfectMatching::BlossomIt "BlossomIt" nested class, |
2077 /// which is able to iterate on the nodes of a blossom. |
1961 /// which is able to iterate on the nodes of a blossom. |
2078 /// If the value type is integer, then the dual solution is multiplied |
1962 /// If the value type is integer, then the dual solution is multiplied |
2079 /// by \ref MaxWeightedMatching::dualScale "4". |
1963 /// by \ref MaxWeightedMatching::dualScale "4". |
2080 /// |
1964 /// |
2081 /// \tparam GR The undirected graph type the algorithm runs on. |
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 /// \ref concepts::Graph::EdgeMap "GR::EdgeMap<int>". |
1967 /// \ref concepts::Graph::EdgeMap "GR::EdgeMap<int>". |
2084 #ifdef DOXYGEN |
1968 #ifdef DOXYGEN |
2085 template <typename GR, typename WM> |
1969 template <typename GR, typename WM> |
2086 #else |
1970 #else |
2087 template <typename GR, |
1971 template <typename GR, |
2989 _delta3->prio() : std::numeric_limits<Value>::max(); |
2870 _delta3->prio() : std::numeric_limits<Value>::max(); |
2990 |
2871 |
2991 Value d4 = !_delta4->empty() ? |
2872 Value d4 = !_delta4->empty() ? |
2992 _delta4->prio() : std::numeric_limits<Value>::max(); |
2873 _delta4->prio() : std::numeric_limits<Value>::max(); |
2993 |
2874 |
2994 _delta_sum = d2; OpType ot = D2; |
2875 _delta_sum = d3; OpType ot = D3; |
2995 if (d3 < _delta_sum) { _delta_sum = d3; ot = D3; } |
2876 if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; } |
2996 if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; } |
2877 if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; } |
2997 |
2878 |
2998 if (_delta_sum == std::numeric_limits<Value>::max()) { |
2879 if (_delta_sum == std::numeric_limits<Value>::max()) { |
2999 return false; |
2880 return false; |
3000 } |
2881 } |
3072 for (NodeIt n(_graph); n != INVALID; ++n) { |
2953 for (NodeIt n(_graph); n != INVALID; ++n) { |
3073 if ((*_matching)[n] != INVALID) { |
2954 if ((*_matching)[n] != INVALID) { |
3074 sum += _weight[(*_matching)[n]]; |
2955 sum += _weight[(*_matching)[n]]; |
3075 } |
2956 } |
3076 } |
2957 } |
3077 return sum /= 2; |
2958 return sum / 2; |
3078 } |
2959 } |
3079 |
2960 |
3080 /// \brief Return \c true if the given edge is in the matching. |
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 /// matching. |
2964 /// matching. |
3084 /// |
2965 /// |
3085 /// \pre Either run() or start() must be called before using this function. |
2966 /// \pre Either run() or start() must be called before using this function. |
3086 bool matching(const Edge& edge) const { |
2967 bool matching(const Edge& edge) const { |
3087 return static_cast<const Edge&>((*_matching)[_graph.u(edge)]) == edge; |
2968 return static_cast<const Edge&>((*_matching)[_graph.u(edge)]) == edge; |
3088 } |
2969 } |
3089 |
2970 |
3090 /// \brief Return the matching arc (or edge) incident to the given node. |
2971 /// \brief Return the matching arc (or edge) incident to the given node. |
3091 /// |
2972 /// |
3092 /// This function returns the matching arc (or edge) incident to the |
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 /// not covered by the matching. |
2975 /// not covered by the matching. |
3095 /// |
2976 /// |
3096 /// \pre Either run() or start() must be called before using this function. |
2977 /// \pre Either run() or start() must be called before using this function. |
3097 Arc matching(const Node& node) const { |
2978 Arc matching(const Node& node) const { |
3098 return (*_matching)[node]; |
2979 return (*_matching)[node]; |
3181 return _blossom_potential[k].value; |
3062 return _blossom_potential[k].value; |
3182 } |
3063 } |
3183 |
3064 |
3184 /// \brief Iterator for obtaining the nodes of a blossom. |
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 /// given blossom. It lists a subset of the nodes. |
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 /// MaxWeightedPerfectMatching class and execute it. |
3070 /// MaxWeightedPerfectMatching class and execute it. |
3190 class BlossomIt { |
3071 class BlossomIt { |
3191 public: |
3072 public: |
3192 |
3073 |
3193 /// \brief Constructor. |
3074 /// \brief Constructor. |
3194 /// |
3075 /// |
3195 /// Constructor to get the nodes of the given variable. |
3076 /// Constructor to get the nodes of the given variable. |
3196 /// |
3077 /// |
3197 /// \pre Either \ref MaxWeightedPerfectMatching::run() "algorithm.run()" |
3078 /// \pre Either \ref MaxWeightedPerfectMatching::run() "algorithm.run()" |
3198 /// or \ref MaxWeightedPerfectMatching::start() "algorithm.start()" |
3079 /// or \ref MaxWeightedPerfectMatching::start() "algorithm.start()" |
3199 /// must be called before initializing this iterator. |
3080 /// must be called before initializing this iterator. |
3200 BlossomIt(const MaxWeightedPerfectMatching& algorithm, int variable) |
3081 BlossomIt(const MaxWeightedPerfectMatching& algorithm, int variable) |
3201 : _algorithm(&algorithm) |
3082 : _algorithm(&algorithm) |
3202 { |
3083 { |
3203 _index = _algorithm->_blossom_potential[variable].begin; |
3084 _index = _algorithm->_blossom_potential[variable].begin; |