Ticket #168: bpm_upgrade.patch
File bpm_upgrade.patch, 81.2 KB (added by , 13 years ago) |
---|
-
lemon/hopcroft_karp.h
# HG changeset patch # User Daniel Poroszkai <poroszd@inf.elte.hu> # Date 1330129198 -3600 # Node ID 6cd8dd8bd86a11f12ed8c617e6af6d5256026603 # Parent 4928e7fa0a4062b96c9fada8abd40d3c7041bcae Hopcroft-Karp improved diff --git a/lemon/hopcroft_karp.h b/lemon/hopcroft_karp.h
a b 2 2 * 3 3 * This file is a part of LEMON, a generic C++ optimization library. 4 4 * 5 * Copyright (C) 2003-201 15 * Copyright (C) 2003-2012 6 6 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport 7 7 * (Egervary Research Group on Combinatorial Optimization, EGRES). 8 8 * … … 20 20 #define LEMON_HOPCROFT_KARP_H 21 21 22 22 #include <lemon/core.h> 23 #include <vector> 24 #include <queue> 23 25 #include <list> 24 26 25 27 /// \ingroup matching … … 35 37 template <typename BPG> 36 38 class HopcroftKarp { 37 39 public: 38 40 /// Type of the bipartite graph 39 41 typedef BPG BpGraph; 40 42 /// Type of the matching map 41 43 typedef typename BPG::template NodeMap<typename BPG::Edge> MatchingMap; … … 81 83 82 84 /// \brief Sets the matching map 83 85 /// 84 /// Sets the matching map. 86 /// Sets the matching map. It should contain a valid matching, 87 /// for example the empty matching is fine. 85 88 /// If you don't use this function before calling \ref init(), 86 89 /// an instance will be allocated automatically. 87 90 /// The destructor deallocates this automatically allocated map, 88 91 /// of course. 89 /// This member is not to initialize the algorithm with a valid90 /// matching; use \ref matchingInit() instead.91 92 /// 92 93 /// \return <tt>(*this)</tt> 93 94 HopcroftKarp& matchingMap(MatchingMap& map) { … … 108 109 109 110 /// \brief Initializes the algorithm 110 111 /// 111 /// Allocates the matching map if necessary, and sets 112 /// to empty matching. 112 /// Allocates the matching map if necessary. 113 113 /// 114 114 /// \pre \ref init() is not called. 115 115 void init() { 116 116 createStructures(); 117 if (!_local_matching) {118 for (NodeIt it(_bpg); it!=INVALID; ++it) {119 _matching->set(it, INVALID);120 }121 }122 117 } 123 118 124 119 /// \brief Smarter initialization of the matching. … … 131 126 /// 132 127 /// \note heuristicInit() replaces init() regarding the preconditions. 133 128 int heuristicInit() { 134 init();129 createStructures(); 135 130 int size = 0; 136 131 137 for (BlueNodeIt b_it(_bpg); b_it!=INVALID; ++b_it) { 138 if ((*_matching)[b_it] == INVALID) { 139 bool matched = false; 140 for (IncEdgeIt inc(_bpg, b_it); inc != INVALID && !matched; ++inc) { 141 if ((*_matching)[_bpg.redNode(inc)] == INVALID) { 142 _matching->set(_bpg.redNode(inc), inc); 143 _matching->set(b_it, inc); 144 matched = true; 145 ++size; 146 } 132 for (BlueNodeIt b(_bpg); b!=INVALID; ++b) { 133 bool matched = false; 134 for (OutArcIt oi(_bpg, b); oi != INVALID && !matched; ++oi) { 135 if ((*_matching)[_bpg.target(oi)] == INVALID) { 136 _matching->set(_bpg.source(oi), oi); 137 _matching->set(_bpg.target(oi), oi); 138 matched = true; 139 ++size; 147 140 } 148 141 } 149 142 } 143 150 144 return size; 151 145 } 152 146 … … 157 151 /// in which an edge is true if it is in the matching. 158 152 /// 159 153 /// If the given matching is invalid, some edges will be left out. 160 /// 154 /// 161 155 /// \return \c false if the given matching is invalid. 162 156 /// 163 157 /// \note matchingInit() replaces init() regarding the preconditions. 164 158 bool matchingInit(const BoolEdgeMap& matching) { 165 init();159 createStructures(); 166 160 bool valid = true; 167 161 for(EdgeIt it(_bpg); it!=INVALID; ++it) { 168 162 if (matching[it]) { … … 200 194 /// 201 195 /// \pre \ref init() must be called before using this function. 202 196 int augment() { 203 IntNodeMap level(_bpg, -1); 204 std::list<RedNode> act_rednodes; 205 std::list<BlueNode> act_bluenodes; 197 IntBlueNodeMap level(_bpg, -2); 198 std::vector<RedNode> act_red_nodes; 206 199 207 200 for (RedNodeIt r_it(_bpg); r_it != INVALID; ++r_it) { 208 201 if ((*_matching)[r_it] == INVALID) { 209 act_red nodes.push_front(r_it);202 act_red_nodes.push_back(r_it); 210 203 } 211 204 } 212 213 // Raise this flag when a shortest augmenting path is found. 214 bool path_found = false; 215 // This nodelist will contain the end nodes of the possible 205 // This vector will contain the end nodes of the possible 216 206 // augmenting paths. 217 std:: list<BlueNode> path_ends;207 std::vector<BlueNode> path_ends; 218 208 219 209 // Layer counter 220 int cur_level = 0;210 int cur_level = -1; 221 211 222 212 // Starting from the unmatched red nodes search for unmatched 223 213 // blue nodes, using Bfs (but only on alternating paths). 224 while (!path_found) { 225 while (!act_rednodes.empty()) { 226 RedNode red = act_rednodes.front(); 227 act_rednodes.pop_front(); 228 level[red] = cur_level; 229 for (IncEdgeIt it(_bpg, red); it!=INVALID; ++it) { 230 BlueNode blue(_bpg.blueNode(it)); 231 if (level[blue] == -1) { 232 act_bluenodes.push_front(blue); 233 level[blue] = cur_level + 1; 234 path_found |= ((*_matching)[blue] == INVALID); 214 // Blue nodes with level==-2 are unreached. 215 while (path_ends.empty()) { 216 cur_level += 2; 217 std::vector<RedNode> next_red_nodes; 218 for (typename std::vector<RedNode>::iterator r=act_red_nodes.begin(); 219 r != act_red_nodes.end(); ++r) { 220 for (OutArcIt it(_bpg, *r); it!=INVALID; ++it) { 221 BlueNode blue(_bpg.asBlueNodeUnsafe(_bpg.target(it))); 222 if (level[blue] != -2) continue; 223 level[blue] = cur_level; 224 if ((*_matching)[blue] == INVALID) { 225 path_ends.push_back(blue); 226 } else { 227 next_red_nodes 228 .push_back(_bpg.redNode((*_matching)[blue])); 235 229 } 236 230 } 237 231 } 238 232 239 if (!path_found) { 240 if (act_bluenodes.empty()) return -1; 241 while (!act_bluenodes.empty()) { 242 BlueNode blue = act_bluenodes.front(); 243 act_bluenodes.pop_front(); 244 RedNode red = _bpg.redNode((*_matching)[blue]); 245 act_rednodes.push_front(red); 246 } 247 } else { 248 for (typename std::list<BlueNode>::iterator it=act_bluenodes.begin(); 249 it != act_bluenodes.end(); 250 ++it) { 251 if ((*_matching)[*it] == INVALID) path_ends.push_front(*it); 252 } 253 // Now path_ends contains nodes that are ends of 254 // shortest alternating paths. 255 } 256 cur_level += 2; 233 if (path_ends.empty() && next_red_nodes.empty()) 234 return -1; 235 236 act_red_nodes.swap(next_red_nodes); 237 next_red_nodes.clear(); 257 238 } 258 259 // The current_edge structure assures that edges iterated at most once 260 typename BpGraph::template BlueNodeMap<IncEdgeIt*> current_edge(_bpg, 0); 239 // length of found augmenting paths 240 int max_level = cur_level; 261 241 262 242 // Stack for the Dfs 263 std::list<BlueNode> stack; 264 243 std::vector<OutArcIt> stack; 244 // Raise this flag when a shortest augmenting path is found. 245 bool path_found; 265 246 // Searching backward, starting from the previously found 266 247 // blue nodes, we build vertex disjoint alternating paths. 267 typename std::list<BlueNode>::iterator n_it = path_ends.begin(); 248 // Blue nodes with level!=-2 are unreached. 249 typename std::vector<BlueNode>::iterator end = path_ends.begin(); 268 250 269 while (n_it != path_ends.end()) { 270 stack.push_front(*n_it); 251 while (end != path_ends.end()) { 252 stack.clear(); 253 stack.push_back(OutArcIt(_bpg, *end)); 254 cur_level = max_level; 271 255 path_found = false; 256 272 257 while (!stack.empty() && !path_found) { 273 BlueNode b = stack.front(); 274 if (current_edge[b] == 0) { 275 current_edge[b] = new IncEdgeIt(_bpg, b); 258 OutArcIt &o = stack.back(); 259 // Search an arc that goes one level higher. 260 // If none found, retreat. 261 while (o != INVALID && 262 (*_matching)[_bpg.target(o)] != INVALID && 263 level[_bpg.blueNode((*_matching)[_bpg.target(o)])] != 264 cur_level - 2) { 265 ++o; 266 } 267 if (o==INVALID) { 268 stack.pop_back(); 269 cur_level += 2; 270 } else if ((*_matching)[_bpg.target(o)] == INVALID) { 271 path_found = true; 276 272 } else { 277 ++(*current_edge[b]); 278 } 279 while (*current_edge[b] != INVALID && 280 level[_bpg.redNode(*current_edge[b])] != level[b] - 1) { 281 ++(*current_edge[b]); 282 } 283 if (*current_edge[b] == INVALID) { 284 stack.pop_front(); 285 } else { 286 RedNode r = _bpg.redNode(*current_edge[b]); 287 if ((*_matching)[r] == INVALID) { 288 path_found = true; 289 } else { 290 level[r] = -1; // Do not visit this node again 291 stack.push_front(_bpg.blueNode((*_matching)[r])); 292 } 273 BlueNode b = _bpg.blueNode((*_matching)[_bpg.target(o)]); 274 stack.push_back(OutArcIt(_bpg, b)); 275 level[b] = -2; // Do not visit this node again 276 cur_level -= 2; 293 277 } 294 278 } 295 279 if (path_found) { 296 BlueNode next = *n_it; 297 RedNode r; 298 BlueNode b; 299 Edge e; 300 while (next != INVALID) { 301 e = *current_edge[next]; 302 b = next; 303 r = _bpg.redNode(e); 304 next = (*_matching)[r] == INVALID ? 305 INVALID :_bpg.blueNode((*_matching)[r]); 306 _matching->set(b, e); 307 _matching->set(r, e); 280 OutArcIt o; 281 while (!stack.empty()) { 282 o = stack.back(); 283 _matching->set(_bpg.source(o), o); 284 _matching->set(_bpg.target(o), o); 285 stack.pop_back(); 308 286 } 309 287 } 310 311 stack.clear(); 312 ++n_it; 288 ++end; 313 289 } 314 return cur_level - 1;290 return max_level; 315 291 } 316 292 317 293 /// \brief Executes the algorithm … … 339 315 /// \brief Size of the matching 340 316 /// 341 317 /// Returns the size of the current matching. 318 /// Function runs in \f$O(n)\f$ time. 342 319 int matchingSize() const { 343 320 int size = 0; 344 321 for (RedNodeIt it(_bpg); it!=INVALID; ++it) { … … 372 349 return (*_matching)[n] != INVALID ? 373 350 _bpg.oppositeNode(n, (*_matching)[n]) : INVALID; 374 351 } 352 353 /// \brief Query a vertex cover. 354 /// 355 /// This function finds a vertex cover based on the current matching, 356 /// and in \c cover sets node in the cover \c true and \c false 357 /// otherwise. In case of maximal matching, this will be a minimal 358 /// vertex cover. 359 /// Function runs in \f$O(n + e)\f$ time. 360 /// 361 /// \return The number of nodes in the cover. 362 int cover(BoolNodeMap &cover) const { 363 int count = 0; 364 std::queue<RedNode> q; 365 366 for (BlueNodeIt b(_bpg); b!=INVALID; ++b) 367 cover[b] = false; 368 369 for (RedNodeIt r(_bpg); r!=INVALID; ++r) { 370 if ((*_matching)[r] == INVALID) { 371 cover[r] = false; 372 q.push(r); 373 } else { 374 cover[r] = true; 375 ++count; 376 } 377 } 378 379 while (!q.empty()) { 380 for (OutArcIt it(_bpg, q.front()); it!=INVALID; ++it) { 381 BlueNode blue(_bpg.asBlueNodeUnsafe(_bpg.target(it))); 382 if (cover[blue]) continue; 383 384 cover[blue] = true; 385 ++count; 386 387 if ((*_matching)[blue] != INVALID) { 388 RedNode nr = _bpg.redNode((*_matching)[blue]); 389 q.push(nr); 390 cover[nr] = false; 391 --count; 392 } 393 } 394 q.pop(); 395 } 396 397 return count; 398 } 399 400 /// \brief Query a barrier of red nodes. 401 /// 402 /// This function tries to create a barrier of red nodes, which should: 403 /// - contain every unmatched red nodes, 404 /// - contain exactly that many matched red nodes as much blue nodes 405 /// there are in its neighbor set. 406 /// 407 /// Barrier exist if and only if the matching is maximal. If exist, the 408 /// map is set true for nodes of the barrier, and false for the others. 409 /// 410 /// Function runs in \f$O(n + e)\f$ time. 411 /// 412 /// \return \c true if a barrier found. 413 bool redBarrier(BoolRedNodeMap &barrier) const { 414 bool exist = true; 415 std::queue<RedNode> q; 416 417 for (RedNodeIt r(_bpg); r!=INVALID; ++r) { 418 if ((*_matching)[r] == INVALID) { 419 barrier[r] = true; 420 q.push(r); 421 } else { 422 barrier[r] = false; 423 } 424 } 425 426 while (!q.empty() && exist) { 427 for (OutArcIt it(_bpg, q.front()); it!=INVALID; ++it) { 428 BlueNode blue(_bpg.asBlueNodeUnsafe(_bpg.target(it))); 429 430 if (exist &= ((*_matching)[blue] != INVALID)) { 431 RedNode r = _bpg.redNode((*_matching)[blue]); 432 if (!barrier[r]) { 433 q.push(r); 434 barrier[r] = true; 435 } 436 } 437 } 438 q.pop(); 439 } 440 441 return exist; 442 } 443 444 /// \brief Query a barrier of blue nodes. 445 /// 446 /// This function tries to create a barrier of blue nodes, which should: 447 /// - contain every unmatched blue nodes, 448 /// - contain exactly that many matched blue nodes as much red nodes 449 /// there are in its neighbor set. 450 /// 451 /// Barrier exist if and only if the matching is maximal. If exist, the 452 /// map is set true for nodes of the barrier, and false for the others. 453 /// 454 /// Function runs in \f$O(n + e)\f$ time. 455 /// 456 /// \return \c true if a barrier found. 457 bool blueBarrier(BoolBlueNodeMap &barrier) const { 458 bool exist = true; 459 std::queue<BlueNode> q; 460 461 for (BlueNodeIt b(_bpg); b!=INVALID; ++b) { 462 if ((*_matching)[b] == INVALID) { 463 barrier[b] = true; 464 q.push(b); 465 } else { 466 barrier[b] = false; 467 } 468 } 469 470 while (!q.empty() && exist) { 471 for (OutArcIt it(_bpg, q.front()); it!=INVALID; ++it) { 472 RedNode r(_bpg.asRedNodeUnsafe(_bpg.target(it))); 473 474 if (exist &= ((*_matching)[r] != INVALID)) { 475 BlueNode b = _bpg.blueNode((*_matching)[r]); 476 if (!barrier[b]) { 477 q.push(b); 478 barrier[b] = true; 479 } 480 } 481 } 482 q.pop(); 483 } 484 485 return exist; 486 } 375 487 }; 376 488 377 489 } -
test/hopcroft_karp_test.cc
diff --git a/test/hopcroft_karp_test.cc b/test/hopcroft_karp_test.cc
a b 2 2 * 3 3 * This file is a part of LEMON, a generic C++ optimization library. 4 4 * 5 * Copyright (C) 2003-201 15 * Copyright (C) 2003-2012 6 6 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport 7 7 * (Egervary Research Group on Combinatorial Optimization, EGRES). 8 8 * … … 104 104 BoolEdgeMap init_matching(graph); 105 105 HK hk_test(graph); 106 106 const HK& const_hk_test = hk_test; 107 HK::MatchingMap matching_map(graph); 107 HK::MatchingMap matching_map(graph, INVALID); 108 BoolNodeMap cov(graph); 109 int cover_size; 110 BoolRedNodeMap rb(graph); 111 BoolBlueNodeMap bb(graph); 108 112 109 113 hk_test.matchingMap(matching_map); 110 114 hk_test.init(); … … 123 127 const_hk_test.mate(node); 124 128 const HK::MatchingMap& max_matching = const_hk_test.matchingMap(); 125 129 edge = max_matching[node]; 130 cover_size = const_hk_test.cover(cov); 131 const_hk_test.redBarrier(rb); 132 const_hk_test.blueBarrier(bb); 126 133 } 127 134 128 135 typedef SmartBpGraph BpGraph; … … 144 151 145 152 } 146 153 147 void checkMatchingOptimality(const BpGraph& bpg, const HK& hk) { 148 list<RedNode> reds; 149 list<BlueNode> blues; 150 BoolBlueNodeMap reached(bpg, false); 154 void checkCover(const BpGraph& bpg, const HK& hk) { 155 int cov_size; 156 BoolNodeMap cov(bpg); 157 cov_size = hk.cover(cov); 158 BoolEdgeMap covered(bpg, false); 159 for (NodeIt n(bpg); n!=INVALID; ++n) { 160 for (IncEdgeIt e(bpg, n); e!=INVALID; ++e) { 161 covered[e] = true; 162 } 163 } 164 bool all = true; 165 for (EdgeIt e(bpg); e!=INVALID; ++e) { 166 all&=covered[e]; 167 } 151 168 152 for (RedNodeIt it(bpg); it != INVALID; ++it) { 153 if (hk.matching(it) == INVALID) { 154 reds.push_front(it); 169 check(all, "Invalid cover"); 170 check(hk.matchingSize() == cov_size, "Suboptimal matching."); 171 } 172 173 void checkBarriers(const BpGraph& bpg, const HK& hk) { 174 BoolRedNodeMap rb(bpg); 175 BoolBlueNodeMap bb(bpg); 176 177 check(hk.redBarrier(rb), "Suboptimal matching."); 178 check(hk.blueBarrier(bb), "Suboptimal matching."); 179 180 BoolNodeMap reached(bpg, false); 181 182 for (RedNodeIt r(bpg); r!=INVALID; ++r) { 183 if (!rb[r]) continue; 184 for (OutArcIt a(bpg, r); a!=INVALID; ++a) { 185 reached.set(bpg.target(a), true); 186 } 187 } 188 for (BlueNodeIt b(bpg); b!=INVALID; ++b) { 189 if (!bb[b]) continue; 190 for (OutArcIt a(bpg, b); a!=INVALID; ++a) { 191 reached.set(bpg.target(a), true); 155 192 } 156 193 } 157 194 158 while (!reds.empty()) { 159 while (!reds.empty()) { 160 RedNode red = reds.front(); 161 reds.pop_front(); 162 for (IncEdgeIt it(bpg, red); it != INVALID; ++it) { 163 BlueNode blue = bpg.blueNode(it); 164 if (!reached[blue]) { 165 blues.push_front(blue); 166 reached.set(blue, true); 167 } 168 } 169 } 170 171 while (!blues.empty()) { 172 BlueNode blue = blues.front(); 173 blues.pop_front(); 174 check(hk.matching(blue) != INVALID, "Suboptimal matching"); 175 reds.push_front(bpg.asRedNode(hk.mate(blue))); 176 } 195 int ur = 0, ub = 0, rc = 0, bc = 0; 196 for (RedNodeIt r(bpg); r!=INVALID; ++r) { 197 if (rb[r]) ++rc; 198 if (hk.mate(r) == INVALID) ++ur; 199 if (reached[r]) --bc; 177 200 } 178 201 202 for (BlueNodeIt b(bpg); b!=INVALID; ++b) { 203 if (bb[b]) ++bc; 204 if (hk.mate(b) == INVALID) ++ub; 205 if (reached[b]) --rc; 206 } 207 208 check(rc == ur, "Red barrier is wrong."); 209 check(bc == ub, "Blue barrier is wrong."); 179 210 } 180 211 212 void checkMatching(const BpGraph& bpg, const HK& hk) { 213 checkMatchingValidity(bpg, hk); 214 checkBarriers(bpg, hk); 215 checkCover(bpg, hk); 216 } 217 218 181 219 int main() { 182 220 // Checking hard wired extremal graphs 183 221 for (int i=0; i<lgfn; ++i) { … … 186 224 bpGraphReader(bpg, lgfs).run(); 187 225 HK hk(bpg); 188 226 hk.run(); 189 checkMatchingValidity(bpg, hk); 190 checkMatchingOptimality(bpg, hk); 227 checkMatching(bpg, hk); 191 228 } 192 229 193 230 // Checking some use cases … … 199 236 hk.init(); 200 237 hk.augment(); 201 238 hk.start(); 202 checkMatchingValidity(bpg, hk); 203 checkMatchingOptimality(bpg, hk); 239 checkMatching(bpg, hk); 204 240 } 205 241 { 206 242 HK hk(bpg); 207 243 hk.heuristicInit(); 208 244 hk.augment(); 209 245 hk.start(); 210 checkMatchingValidity(bpg, hk); 211 checkMatchingOptimality(bpg, hk); 246 checkMatching(bpg, hk); 212 247 } 213 248 { 214 249 HK hk(bpg); 215 HK::MatchingMap m(bpg );250 HK::MatchingMap m(bpg, INVALID); 216 251 hk.matchingMap(m); 217 252 hk.run(); 218 checkMatchingValidity(bpg, hk); 219 checkMatchingOptimality(bpg, hk); 253 checkMatching(bpg, hk); 220 254 } 221 255 { 222 256 HK hk(bpg); … … 226 260 init_matching[e] = true; 227 261 check(hk.matchingInit(init_matching), "Wrong result from matchingInit()"); 228 262 hk.start(); 229 checkMatchingValidity(bpg, hk); 230 checkMatchingOptimality(bpg, hk); 263 checkMatching(bpg, hk); 231 264 } 232 265 { 233 266 HK hk(bpg); 234 HK::MatchingMap m(bpg );267 HK::MatchingMap m(bpg, INVALID); 235 268 hk.matchingMap(m); 236 269 BpGraph::EdgeMap<bool> init_matching(bpg, true); 237 270 check(!hk.matchingInit(init_matching), "Wrong result from matchingInit()"); 238 271 hk.start(); 239 checkMatchingValidity(bpg, hk); 240 checkMatchingOptimality(bpg, hk); 272 checkMatching(bpg, hk); 241 273 } 242 243 274 // Checking some random generated graphs 244 275 const int random_test_n = 10; 245 const int max_red_n = 100 ;276 const int max_red_n = 1000; 246 277 const int min_red_n = 10; 247 const int max_blue_n = 100 ;278 const int max_blue_n = 1000; 248 279 const int min_blue_n = 10; 249 280 for (int i=0; i<random_test_n; ++i) { 250 281 int red_n = rnd.integer(min_red_n, max_red_n); … … 259 290 } 260 291 for (SmartBpGraph::RedNodeIt r(bpg); r!=INVALID; ++r) { 261 292 for (SmartBpGraph::BlueNodeIt b(bpg); b!=INVALID; ++b) { 262 if (rnd() < prob ) {293 if (rnd() < prob/10.) { 263 294 bpg.addEdge(r,b); 264 295 } 265 296 } … … 267 298 268 299 HK hk(bpg); 269 300 hk.run(); 270 checkMatchingValidity(bpg, hk); 271 checkMatchingOptimality(bpg, hk); 301 checkMatching(bpg, hk); 272 302 } 273 303 274 304 -
lemon/elevator.h
# HG changeset patch # User Daniel Poroszkai <poroszd@inf.elte.hu> # Date 1329865937 -3600 # Node ID b5f7d8e60bc7f599bf0e2bf60bda1e05cf3c6d47 # Parent dbf9516f49c509e59a9c73da18eb868cfafb751c Fix initialization bug in elevator diff --git a/lemon/elevator.h b/lemon/elevator.h
a b 2 2 * 3 3 * This file is a part of LEMON, a generic C++ optimization library. 4 4 * 5 * Copyright (C) 2003-20 095 * Copyright (C) 2003-2012 6 6 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport 7 7 * (Egervary Research Group on Combinatorial Optimization, EGRES). 8 8 * … … 108 108 Elevator(const GR &graph,int max_level) : 109 109 _g(graph), 110 110 _max_level(max_level), 111 _item_num( _max_level),111 _item_num(countItems<GR, Item>(graph)), 112 112 _where(graph), 113 113 _level(graph,0), 114 _items(_ max_level),114 _items(_item_num), 115 115 _first(_max_level+2), 116 116 _last_active(_max_level+2), 117 117 _highest_active(-1) {} … … 526 526 ///\param max_level The maximum allowed level. 527 527 ///Set the range of the possible labels to <tt>[0..max_level]</tt>. 528 528 LinkedElevator(const GR& graph, int max_level) 529 : _graph(graph), _max_level(max_level), _item_num(_max_level), 529 : _graph(graph), _max_level(max_level), 530 _item_num(countItems<GR, Item>(graph)), 530 531 _first(_max_level + 1), _last(_max_level + 1), 531 532 _prev(graph), _next(graph), 532 533 _highest_active(-1), _level(graph), _active(graph) {} -
new file lemon/preflow_bp_matching.h
# HG changeset patch # User Daniel Poroszkai <poroszd@inf.elte.hu> # Date 1330129742 -3600 # Node ID 29f2b9488e4ea58b107564fa272ae38b402b2c1e # Parent b5f7d8e60bc7f599bf0e2bf60bda1e05cf3c6d47 Preflow based bipartite matching diff --git a/lemon/preflow_bp_matching.h b/lemon/preflow_bp_matching.h new file mode 100644
- + 1 /* -*- mode: C++; indent-tabs-mode: nil; -*- 2 * 3 * This file is a part of LEMON, a generic C++ optimization library. 4 * 5 * Copyright (C) 2003-2012 6 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport 7 * (Egervary Research Group on Combinatorial Optimization, EGRES). 8 * 9 * Permission to use, modify and distribute this software is granted 10 * provided that this copyright notice appears in all copies. For 11 * precise terms see the accompanying LICENSE file. 12 * 13 * This software is provided "AS IS" with no warranty of any kind, 14 * express or implied, and with no claim as to its suitability for any 15 * purpose. 16 * 17 */ 18 19 #ifndef LEMON_PREFLOW_BP_MATCHING_H 20 #define LEMON_PREFLOW_BP_MATCHING_H 21 22 #include <lemon/core.h> 23 #include <lemon/elevator.h> 24 #include <vector> 25 #include <list> 26 #include <queue> 27 28 /// \ingroup matching 29 /// \file 30 /// \brief A preflow based bipartite matching algorithm. 31 namespace lemon { 32 33 /// \brief A preflow based bipartite matching algorithm 34 /// 35 /// Finds maximum cardinality matching in a given bipartite 36 /// graph using specialized preflow algorithm. 37 template <typename BPG> 38 class PreflowBpMatching { 39 public: 40 /// Type of the bipartite graph 41 typedef BPG BpGraph; 42 /// Type of the matching map 43 typedef typename BPG::template NodeMap<typename BPG::Edge> MatchingMap; 44 /// Type of the elevator. 45 typedef Elevator<BPG, typename BPG::BlueNode> Elevator_t; 46 47 private: 48 TEMPLATE_BPGRAPH_TYPEDEFS(BPG); 49 50 protected: 51 typedef typename BPG::template RedNodeMap<typename BPG::Arc> Preflow; 52 typedef std::list<typename BPG::Arc> CurrentArcSet; 53 typedef typename 54 std::list<typename BPG::Arc>::iterator CurrentArcIterator; 55 56 const BpGraph& _bpg; 57 MatchingMap* _matching; 58 bool _local_matching; 59 Elevator_t* _elevator; 60 bool _local_elevator; 61 typename BPG::template BlueNodeMap<CurrentArcSet> _cas; 62 typename BPG::template BlueNodeMap<CurrentArcIterator> _ca; 63 Preflow _flow; 64 65 void createStructures() { 66 if (!_matching) { 67 _matching = new MatchingMap(_bpg, INVALID); 68 _local_matching = true; 69 } 70 if (!_elevator) { 71 _elevator = new Elevator_t(_bpg, std::min(countRedNodes(_bpg) + 1, 72 countBlueNodes(_bpg))); 73 _local_elevator = true; 74 } 75 } 76 77 void destroyStructures() { 78 if (_local_matching) { 79 delete _matching; 80 } 81 if (_local_elevator) { 82 delete _elevator; 83 } 84 } 85 86 /// \brief Initialize the elevator 87 /// 88 /// When this function is called, the elevator is initialized with 89 /// the exact distance labels corresponding the current matching. 90 void initElevator() { 91 92 // Set the initial flow: 93 // choose arbitrary outgoing arc for every unmatched red node. 94 for (RedNodeIt r(_bpg); r!=INVALID; ++r) { 95 Arc a = (*_matching)[r] != INVALID ? 96 _bpg.direct((*_matching)[r], true) : OutArcIt(_bpg, r); 97 if (a != INVALID) { 98 _flow[r] = a; 99 _cas[_bpg.asBlueNodeUnsafe(_bpg.target(a))].push_front(a); 100 } 101 } 102 // Now we start a backward Bfs from the unmatched blue nodes in the 103 // residual graph, to supply the exact distance labels to the elevator. 104 std::vector<BlueNode> act_blue_nodes, next_blue_nodes; 105 BoolBlueNodeMap reached(_bpg, false); 106 107 for (BlueNodeIt b(_bpg); b != INVALID; ++b) { 108 if (_cas[b].empty()) { 109 act_blue_nodes.push_back(b); 110 reached.set(b, true); 111 } 112 } 113 114 _elevator->initStart(); 115 while (!act_blue_nodes.empty()) { 116 for (typename std::vector<BlueNode>::iterator b_it = 117 act_blue_nodes.begin(); 118 b_it != act_blue_nodes.end(); ++b_it) { 119 120 _elevator->initAddItem(*b_it); 121 for (OutArcIt oi(_bpg, *b_it); oi!=INVALID; ++oi) { 122 RedNode r = _bpg.asRedNodeUnsafe(_bpg.target(oi)); 123 BlueNode b = _bpg.asBlueNodeUnsafe(_bpg.target(_flow[r])); 124 if (b != *b_it && !reached[b]) { 125 next_blue_nodes.push_back(b); 126 reached.set(b, true); 127 } 128 } 129 } 130 _elevator->initNewLevel(); 131 next_blue_nodes.swap(act_blue_nodes); 132 next_blue_nodes.clear(); 133 } 134 _elevator->initFinish(); 135 136 // Activate blue nodes with excess, i. e. those 137 // which get flow from more than one red node. 138 for (BlueNodeIt b(_bpg); b != INVALID; ++b) { 139 if (_cas[b].size() > 1 && (*_elevator)[b] < _elevator->maxLevel()) { 140 _elevator->activate(b); 141 } 142 143 _ca[b] = _cas[b].begin(); 144 } 145 } 146 147 PreflowBpMatching() {} 148 149 public: 150 151 /// \brief Constructor 152 /// 153 /// Constructs the class on the given bipartite graph. 154 PreflowBpMatching(const BpGraph& bpg) : _bpg(bpg), 155 _matching(0), 156 _local_matching(false), 157 _elevator(0), 158 _local_elevator(false), 159 _cas(bpg), 160 _ca(bpg), 161 _flow(bpg, INVALID) 162 {} 163 164 /// \brief Destructor 165 /// 166 /// Destructor. 167 ~PreflowBpMatching() { 168 destroyStructures(); 169 } 170 171 /// \brief Sets the matching map 172 /// 173 /// Sets the matching map. It should contain a valid matching, 174 /// for example the empty matching is fine. 175 /// If you don't use this function before calling \ref init(), 176 /// an instance will be allocated automatically. 177 /// The destructor deallocates this automatically allocated map, 178 /// of course. 179 /// 180 /// \return <tt>(*this)</tt> 181 PreflowBpMatching& matchingMap(MatchingMap& map) { 182 _matching = ↦ 183 _local_matching = false; 184 return *this; 185 } 186 187 /// \brief Returns a const reference to the matching map. 188 /// 189 /// Returns a const reference to the matching map, which contains 190 /// the matching edge for each node (or INVALID for unmatched nodes). 191 const MatchingMap& matchingMap() const { 192 return *_matching; 193 } 194 195 /// \brief Sets the elevator used by algorithm. 196 /// 197 /// Sets the elevator. 198 /// If you don't use this function before calling \ref init(), 199 /// an instance will be allocated automatically. 200 /// The destructor deallocates this automatically allocated elevator, 201 /// of course. 202 /// \return <tt>(*this)</tt> 203 PreflowBpMatching& elevator(Elevator_t& elevator) { 204 _elevator = &elevator; 205 _local_elevator = false; 206 return *this; 207 } 208 209 /// \brief Returns a const reference to the elevator. 210 /// 211 /// Returns a const reference to the elevator, which keeps track 212 /// of the distance labels of nodes. 213 const Elevator_t& elevator() const { 214 return *_elevator; 215 } 216 217 /// \brief Initializes the algorithm 218 /// 219 /// Allocates the elevator and the matching map if necessary, 220 /// and initializes the elevator. 221 /// 222 /// \pre \ref init() is not called. 223 void init() { 224 createStructures(); 225 initElevator(); 226 } 227 228 /// \brief Smarter initialization of the matching. 229 /// 230 /// Allocates the matching map if necessary, and initializes 231 /// the algorithm with a trivially found matching: iterating 232 /// on every edge, take it in if both endnodes are unmatched. 233 /// 234 /// \return Size of the so found matching. 235 /// 236 /// \note heuristicInit() replaces init() regarding the preconditions. 237 int heuristicInit() { 238 createStructures(); 239 int size = 0; 240 241 for (BlueNodeIt b(_bpg); b!=INVALID; ++b) { 242 bool matched = false; 243 for (OutArcIt oi(_bpg, b); oi != INVALID && !matched; ++oi) { 244 if ((*_matching)[_bpg.target(oi)] == INVALID) { 245 _matching->set(_bpg.source(oi), oi); 246 _matching->set(_bpg.target(oi), oi); 247 matched = true; 248 ++size; 249 } 250 } 251 } 252 253 initElevator(); 254 return size; 255 } 256 257 /// \brief Initialize the matching from a map. 258 /// 259 /// Allocates the matching map if necessary, and 260 /// initializes the algorithm with a matching given as a bool edgemap, 261 /// in which an edge is true if it is in the matching. 262 /// Also initializes the elevator. 263 /// 264 /// \return \c true if the matching is valid. 265 /// 266 /// \note matchingInit() replaces init() regarding the preconditions. 267 bool matchingInit(const BoolEdgeMap& matching) { 268 createStructures(); 269 270 bool valid = true; 271 for(EdgeIt it(_bpg); it!=INVALID; ++it) { 272 if (matching[it]) { 273 Node red = _bpg.redNode(it); 274 Node blue = _bpg.blueNode(it); 275 if ((*_matching)[red] != INVALID || (*_matching)[blue] != INVALID) { 276 valid = false; 277 } else { 278 _matching->set(red, it); 279 _matching->set(blue, it); 280 } 281 } 282 } 283 284 initElevator(); 285 return valid; 286 } 287 288 289 /// \brief Perform a push/relabel operation on a node 290 /// 291 /// Perform a push/relabel operation on a node. 292 /// 293 /// \pre \ref init() is called, and the node is active. 294 void push_relabel(const BlueNode& chosen) { 295 Arc new_flow = INVALID; 296 while (_ca[chosen] != _cas[chosen].end() && new_flow == INVALID) { 297 OutArcIt r_oi(_bpg, _bpg.source(*_ca[chosen])); 298 while (r_oi != INVALID && new_flow == INVALID) { 299 if ((*_elevator)[_bpg.asBlueNodeUnsafe(_bpg.target(r_oi))] == 300 (*_elevator)[chosen] - 1) { 301 new_flow = r_oi; 302 } 303 ++r_oi; 304 } 305 if (new_flow == INVALID) ++_ca[chosen]; 306 } 307 308 // If an incident blue node is found on one level lower, 309 // then perform a push, i. e. flip the flow from the 310 // current arc to 'new_flow'. 311 if (new_flow != INVALID) { 312 _ca[chosen] = _cas[chosen].erase(_ca[chosen]); 313 if (_cas[chosen].size() == 1) { 314 _elevator->deactivate(chosen); 315 } 316 _flow[_bpg.asRedNodeUnsafe(_bpg.source(new_flow))] = new_flow; 317 BlueNode b = _bpg.asBlueNodeUnsafe(_bpg.target(new_flow)); 318 _cas[b].push_front(new_flow); 319 if (_cas[b].size() == 2) { 320 _elevator->activate(b); 321 } 322 } 323 324 // If none found, then relabel. 325 else { 326 int min = _elevator->maxLevel() + 1; 327 for (CurrentArcIterator it1 = _cas[chosen].begin(); 328 it1 != _cas[chosen].end(); 329 ++it1) { 330 for (OutArcIt it2(_bpg, _bpg.asRedNodeUnsafe(_bpg.source(*it1))); 331 it2 != INVALID; 332 ++it2) { 333 BlueNode b = _bpg.asBlueNodeUnsafe(_bpg.target(it2)); 334 if (b != chosen && min > (*_elevator)[b]) { 335 min = (*_elevator)[b]; 336 } 337 } 338 } 339 340 int chosen_level = (*_elevator)[chosen]; 341 if (min >= _elevator->maxLevel() - 1) { 342 _elevator->lift(chosen, _elevator->maxLevel()); 343 _elevator->deactivate(chosen); 344 } else { 345 _elevator->lift(chosen, min+1); 346 } 347 if (_elevator->emptyLevel(chosen_level)) { 348 _elevator->liftToTop(chosen_level); 349 } 350 351 // Restore the current arc iterator. 352 _ca[chosen] = _cas[chosen].begin(); 353 } 354 } 355 356 /// \brief Executes the algorithm 357 /// 358 /// Executes push/relabel operation on the highest active node 359 /// until there is. After start(), writeMatching() should be called to 360 /// remove unnecessary flow and choose a valid maximum matching. 361 /// 362 /// \pre \ref init() must be called before using this function. 363 void start() { 364 while (_elevator->highestActive() != INVALID) 365 push_relabel(_elevator->highestActive()); 366 } 367 368 /// \brief Deduces and writes a matching from the preflow 369 /// 370 /// The underlying preflow algorithm does not use the matching map, 371 /// so before using the matching query functions ( \ref mate(), \ref 372 /// matchingSize(), etc.), you should explicitly call this function 373 /// to write a matching based on the current preflow to the matching 374 /// map. When the preflow results an ambiguous matching (a node has 375 /// incoming flow on more than one arc), it chooses arbitrary one. 376 void writeMatching() { 377 for (NodeIt it(_bpg); it != INVALID; ++it) { 378 _matching->set(it, INVALID); 379 } 380 for (BlueNodeIt b(_bpg); b != INVALID; ++b) { 381 if (!_cas[b].empty()) { 382 Arc a = _cas[b].front(); 383 _matching->set(_bpg.source(a), a); 384 _matching->set(_bpg.target(a), a); 385 } 386 } 387 } 388 389 /// \brief Runs the algorithm. 390 /// 391 /// bp1.run() is just a shorthand for: 392 /// 393 ///\code 394 /// bp1.heuristicInit(); 395 /// bp1.start(); 396 /// bp1.writeMatching(); 397 ///\endcode 398 void run() { 399 heuristicInit(); 400 start(); 401 writeMatching(); 402 } 403 404 /// \brief Size of the matching 405 /// 406 /// Returns the size of the current matching. 407 /// Function runs in \f$O(n)\f$ time. 408 int matchingSize() const { 409 int size = 0; 410 for (RedNodeIt it(_bpg); it!=INVALID; ++it) { 411 if ((*_matching)[it] != INVALID) ++size; 412 } 413 return size; 414 } 415 416 /// \brief Return \c true if the given edge is in the matching. 417 /// 418 /// This function returns \c true if the given edge is in the current 419 /// matching. 420 bool matching(const Edge& edge) const { 421 return edge == (*_matching)[_bpg.redNode(edge)]; 422 } 423 424 /// \brief Return the matching edge incident to the given node. 425 /// 426 /// This function returns the matching edge incident to the 427 /// given node in the current matching or \c INVALID if the node is 428 /// not covered by the matching. 429 Edge matching(const Node& n) const { 430 return (*_matching)[n]; 431 } 432 433 /// \brief Return the mate of the given node. 434 /// 435 /// This function returns the mate of the given node in the current 436 /// matching or \c INVALID if the node is not covered by the matching. 437 Node mate(const Node& n) const { 438 return (*_matching)[n] != INVALID ? 439 _bpg.oppositeNode(n, (*_matching)[n]) : INVALID; 440 } 441 442 /// \brief Query a vertex cover. 443 /// 444 /// This function finds a vertex cover based on the current matching, 445 /// and in \c cover sets node in the cover \c true and \c false 446 /// otherwise. In case of maximal matching, this will be a minimal 447 /// vertex cover. 448 /// Function runs in \f$O(n + e)\f$ time. 449 /// 450 /// \return The number of nodes in the cover. 451 int cover(BoolNodeMap &cover) const { 452 int count = 0; 453 std::queue<RedNode> q; 454 455 for (BlueNodeIt b(_bpg); b!=INVALID; ++b) 456 cover[b] = false; 457 458 for (RedNodeIt r(_bpg); r!=INVALID; ++r) { 459 if ((*_matching)[r] == INVALID) { 460 cover[r] = false; 461 q.push(r); 462 } else { 463 cover[r] = true; 464 ++count; 465 } 466 } 467 468 while (!q.empty()) { 469 for (OutArcIt it(_bpg, q.front()); it!=INVALID; ++it) { 470 BlueNode blue(_bpg.asBlueNodeUnsafe(_bpg.target(it))); 471 if (cover[blue]) continue; 472 473 cover[blue] = true; 474 ++count; 475 476 if ((*_matching)[blue] != INVALID) { 477 RedNode nr = _bpg.redNode((*_matching)[blue]); 478 q.push(nr); 479 cover[nr] = false; 480 --count; 481 } 482 } 483 q.pop(); 484 } 485 486 return count; 487 } 488 489 /// \brief Query a barrier of red nodes. 490 /// 491 /// This function tries to create a barrier of red nodes, which should: 492 /// - contain every unmatched red nodes, 493 /// - contain exactly that many matched red nodes as much blue nodes 494 /// there are in its neighbor set. 495 /// 496 /// Barrier exist if and only if the matching is maximal. If exist, the 497 /// map is set true for nodes of the barrier, and false for the others. 498 /// 499 /// Function runs in \f$O(n + e)\f$ time. 500 /// 501 /// \return \c true if a barrier found. 502 bool redBarrier(BoolRedNodeMap &barrier) const { 503 bool exist = true; 504 std::queue<RedNode> q; 505 506 for (RedNodeIt r(_bpg); r!=INVALID; ++r) { 507 if ((*_matching)[r] == INVALID) { 508 barrier[r] = true; 509 q.push(r); 510 } else { 511 barrier[r] = false; 512 } 513 } 514 515 while (!q.empty() && exist) { 516 for (OutArcIt it(_bpg, q.front()); it!=INVALID; ++it) { 517 BlueNode blue(_bpg.asBlueNodeUnsafe(_bpg.target(it))); 518 519 if (exist &= ((*_matching)[blue] != INVALID)) { 520 RedNode r = _bpg.redNode((*_matching)[blue]); 521 if (!barrier[r]) { 522 q.push(r); 523 barrier[r] = true; 524 } 525 } 526 } 527 q.pop(); 528 } 529 530 return exist; 531 } 532 533 /// \brief Query a barrier of blue nodes. 534 /// 535 /// This function tries to create a barrier of blue nodes, which should: 536 /// - contain every unmatched blue nodes, 537 /// - contain exactly that many matched blue nodes as much red nodes 538 /// there are in its neighbor set. 539 /// 540 /// Barrier exist if and only if the matching is maximal. If exist, the 541 /// map is set true for nodes of the barrier, and false for the others. 542 /// 543 /// Function runs in \f$O(n + e)\f$ time. 544 /// 545 /// \return \c true if a barrier found. 546 bool blueBarrier(BoolBlueNodeMap &barrier) const { 547 bool exist = true; 548 std::queue<BlueNode> q; 549 550 for (BlueNodeIt b(_bpg); b!=INVALID; ++b) { 551 if ((*_matching)[b] == INVALID) { 552 barrier[b] = true; 553 q.push(b); 554 } else { 555 barrier[b] = false; 556 } 557 } 558 559 while (!q.empty() && exist) { 560 for (OutArcIt it(_bpg, q.front()); it!=INVALID; ++it) { 561 RedNode r(_bpg.asRedNodeUnsafe(_bpg.target(it))); 562 563 if (exist &= ((*_matching)[r] != INVALID)) { 564 BlueNode b = _bpg.blueNode((*_matching)[r]); 565 if (!barrier[b]) { 566 q.push(b); 567 barrier[b] = true; 568 } 569 } 570 } 571 q.pop(); 572 } 573 574 return exist; 575 } 576 }; 577 578 } 579 #endif -
test/CMakeLists.txt
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
a b 48 48 path_test 49 49 planarity_test 50 50 preflow_test 51 preflow_bp_matching_test 51 52 radix_sort_test 52 53 random_test 53 54 suurballe_test -
new file test/preflow_bp_matching_test.cc
diff --git a/test/preflow_bp_matching_test.cc b/test/preflow_bp_matching_test.cc new file mode 100644
- + 1 /* -*- mode: C++; indent-tabs-mode: nil; -*- 2 * 3 * This file is a part of LEMON, a generic C++ optimization library. 4 * 5 * Copyright (C) 2003-2012 6 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport 7 * (Egervary Research Group on Combinatorial Optimization, EGRES). 8 * 9 * Permission to use, modify and distribute this software is granted 10 * provided that this copyright notice appears in all copies. For 11 * precise terms see the accompanying LICENSE file. 12 * 13 * This software is provided "AS IS" with no warranty of any kind, 14 * express or implied, and with no claim as to its suitability for any 15 * purpose. 16 * 17 */ 18 19 #include <iostream> 20 #include <sstream> 21 #include <list> 22 23 #include <lemon/random.h> 24 #include <lemon/preflow_bp_matching.h> 25 #include <lemon/smart_graph.h> 26 #include <lemon/concepts/bpgraph.h> 27 #include <lemon/concepts/maps.h> 28 #include <lemon/lgf_reader.h> 29 30 #include "test_tools.h" 31 32 using namespace std; 33 using namespace lemon; 34 35 const int lgfn = 3; 36 const string lgf[lgfn] = { 37 // Empty graph 38 "@red_nodes\n" 39 "label\n" 40 "\n" 41 "@blue_nodes\n" 42 "label\n" 43 "\n" 44 "@edges\n" 45 " label\n" 46 "\n", 47 48 // No red nodes 49 "@red_nodes\n" 50 "label\n" 51 "\n" 52 "@blue_nodes\n" 53 "label\n" 54 "1\n" 55 "@edges\n" 56 " label\n" 57 "\n", 58 59 // No blue nodes 60 "@red_nodes\n" 61 "label\n" 62 "1\n" 63 "@blue_nodes\n" 64 "label\n" 65 "\n" 66 "@edges\n" 67 " label\n" 68 "\n", 69 }; 70 71 const string u_lgf = 72 "@red_nodes\n" 73 "label\n" 74 "1\n" 75 "2\n" 76 "3\n" 77 "@blue_nodes\n" 78 "label\n" 79 "4\n" 80 "5\n" 81 "6\n" 82 "7\n" 83 "@edges\n" 84 "label\n" 85 "1 4 1\n" 86 "1 5 2\n" 87 "2 5 3\n" 88 "2 6 4\n" 89 "2 7 5\n" 90 "3 7 6\n" 91 ; 92 93 94 95 void checkPreflowBpMatchingCompile() { 96 typedef concepts::BpGraph BpGraph; 97 BPGRAPH_TYPEDEFS(BpGraph); 98 typedef PreflowBpMatching<BpGraph> PBM; 99 100 BpGraph graph; 101 RedNode red; 102 BlueNode blue; 103 Node node; 104 Edge edge; 105 BoolEdgeMap init_matching(graph); 106 PBM pbm_test(graph); 107 const PBM& const_pbm_test = pbm_test; 108 PBM::MatchingMap matching_map(graph, INVALID); 109 PBM::Elevator_t elevator(graph); 110 BpGraph::NodeMap<bool> cov(graph); 111 int cover_size; 112 BoolRedNodeMap rb(graph); 113 BoolBlueNodeMap bb(graph); 114 115 pbm_test.matchingMap(matching_map); 116 pbm_test.elevator(elevator); 117 pbm_test.init(); 118 pbm_test.heuristicInit(); 119 pbm_test.matchingInit(init_matching); 120 121 pbm_test.start(); 122 pbm_test.push_relabel(blue); 123 pbm_test.run(); 124 pbm_test.writeMatching(); 125 126 const_pbm_test.matchingSize(); 127 const_pbm_test.matching(node); 128 const_pbm_test.matching(edge); 129 const_pbm_test.mate(red); 130 const_pbm_test.mate(blue); 131 const_pbm_test.mate(node); 132 const PBM::MatchingMap& max_matching = const_pbm_test.matchingMap(); 133 edge = max_matching[node]; 134 const PBM::Elevator_t& elev = const_pbm_test.elevator(); 135 blue = elev.highestActive(); 136 cover_size = const_pbm_test.cover(cov); 137 const_pbm_test.redBarrier(rb); 138 const_pbm_test.blueBarrier(bb); 139 } 140 141 typedef SmartBpGraph BpGraph; 142 typedef PreflowBpMatching<BpGraph> PBM; 143 BPGRAPH_TYPEDEFS(BpGraph); 144 145 void checkMatchingValidity(const BpGraph& bpg, const PBM& pbm) { 146 BoolBlueNodeMap matched(bpg, false); 147 for (RedNodeIt it(bpg); it != INVALID; ++it) { 148 if (pbm.matching(it) != INVALID) { 149 check(pbm.mate(pbm.mate(it)) == it, "Wrong matching"); 150 matched.set(bpg.asBlueNode(pbm.mate(it)), true); 151 } 152 } 153 for (BlueNodeIt it(bpg); it != INVALID; ++it) { 154 // != is used as logical xor here 155 check(matched[it] != (pbm.matching(it) == INVALID), "Wrong matching"); 156 } 157 } 158 159 void checkCover(const BpGraph& bpg, const PBM& pbm) { 160 int cov_size; 161 BoolNodeMap cov(bpg); 162 cov_size = pbm.cover(cov); 163 BoolEdgeMap covered(bpg, false); 164 for (NodeIt n(bpg); n!=INVALID; ++n) { 165 for (IncEdgeIt e(bpg, n); e!=INVALID; ++e) { 166 covered[e] = true; 167 } 168 } 169 bool all = true; 170 for (EdgeIt e(bpg); e!=INVALID; ++e) { 171 all&=covered[e]; 172 } 173 174 check(all, "Invalid cover"); 175 check(pbm.matchingSize() == cov_size, "Suboptimal matching."); 176 } 177 178 void checkBarriers(const BpGraph& bpg, const PBM& pbm) { 179 BoolRedNodeMap rb(bpg); 180 BoolBlueNodeMap bb(bpg); 181 182 check(pbm.redBarrier(rb), "Suboptimal matching."); 183 check(pbm.blueBarrier(bb), "Suboptimal matching."); 184 185 BoolNodeMap reached(bpg, false); 186 187 for (RedNodeIt r(bpg); r!=INVALID; ++r) { 188 if (!rb[r]) continue; 189 for (OutArcIt a(bpg, r); a!=INVALID; ++a) { 190 reached.set(bpg.target(a), true); 191 } 192 } 193 for (BlueNodeIt b(bpg); b!=INVALID; ++b) { 194 if (!bb[b]) continue; 195 for (OutArcIt a(bpg, b); a!=INVALID; ++a) { 196 reached.set(bpg.target(a), true); 197 } 198 } 199 200 int ur = 0, ub = 0, rc = 0, bc = 0; 201 for (RedNodeIt r(bpg); r!=INVALID; ++r) { 202 if (rb[r]) ++rc; 203 if (pbm.mate(r) == INVALID) ++ur; 204 if (reached[r]) --bc; 205 } 206 207 for (BlueNodeIt b(bpg); b!=INVALID; ++b) { 208 if (bb[b]) ++bc; 209 if (pbm.mate(b) == INVALID) ++ub; 210 if (reached[b]) --rc; 211 } 212 213 check(rc == ur, "Red barrier is wrong."); 214 check(bc == ub, "Blue barrier is wrong."); 215 } 216 217 void checkMatching(const BpGraph& bpg, const PBM& pbm) { 218 checkMatchingValidity(bpg, pbm); 219 checkBarriers(bpg, pbm); 220 checkCover(bpg, pbm); 221 } 222 223 224 int main() { 225 // Checking hard wired extremal graphs 226 for (int i=0; i<lgfn; ++i) { 227 BpGraph bpg; 228 istringstream lgfs(lgf[i]); 229 bpGraphReader(bpg, lgfs).run(); 230 PBM pbm(bpg); 231 pbm.run(); 232 checkMatching(bpg, pbm); 233 } 234 235 // Checking some use cases 236 BpGraph bpg; 237 istringstream u_lgfs(u_lgf); 238 bpGraphReader(bpg, u_lgfs).run(); 239 240 { 241 PBM pbm(bpg); 242 pbm.init(); 243 BpGraph::BlueNode b = pbm.elevator().highestActive(); 244 if (b != INVALID) { 245 pbm.push_relabel(b); 246 } 247 pbm.start(); 248 pbm.writeMatching(); 249 check(pbm.elevator().highestActive() == INVALID, "Error in elevator"); 250 checkMatching(bpg, pbm); 251 } 252 253 { 254 PBM pbm(bpg); 255 pbm.heuristicInit(); 256 BpGraph::BlueNode b = pbm.elevator().highestActive(); 257 if (b != INVALID) { 258 pbm.push_relabel(b); 259 } 260 pbm.start(); 261 checkMatching(bpg, pbm); 262 } 263 { 264 PBM pbm(bpg); 265 PBM::MatchingMap m(bpg, INVALID); 266 pbm.matchingMap(m); 267 pbm.run(); 268 checkMatching(bpg, pbm); 269 } 270 { 271 PBM pbm(bpg); 272 BpGraph::EdgeMap<bool> init_matching(bpg, false); 273 BpGraph::Edge e; 274 bpg.first(e); 275 init_matching[e] = true; 276 check(pbm.matchingInit(init_matching), 277 "Wrong result from matchingInit()"); 278 pbm.start(); 279 pbm.writeMatching(); 280 checkMatching(bpg, pbm); 281 } 282 { 283 PBM pbm(bpg); 284 PBM::MatchingMap m(bpg, INVALID); 285 pbm.matchingMap(m); 286 BpGraph::EdgeMap<bool> init_matching(bpg, true); 287 check(!pbm.matchingInit(init_matching), 288 "Wrong result from matchingInit()"); 289 pbm.start(); 290 pbm.writeMatching(); 291 checkMatching(bpg, pbm); 292 } 293 294 // Checking some random generated graphs 295 const int random_test_n = 10; 296 const int max_red_n = 1000; 297 const int min_red_n = 10; 298 const int max_blue_n = 1000; 299 const int min_blue_n = 10; 300 for (int i=0; i<random_test_n; ++i) { 301 int red_n = rnd.integer(min_red_n, max_red_n); 302 int blue_n = rnd.integer(min_blue_n, max_blue_n); 303 double prob = rnd(); 304 305 SmartBpGraph bpg; 306 for (int j=0; j<red_n; ++j) { 307 bpg.addRedNode(); 308 } 309 for (int j=0; j<blue_n; ++j) { 310 bpg.addBlueNode(); 311 } 312 for (SmartBpGraph::RedNodeIt r(bpg); r!=INVALID; ++r) { 313 for (SmartBpGraph::BlueNodeIt b(bpg); b!=INVALID; ++b) { 314 if (rnd() < prob/10.) { 315 bpg.addEdge(r,b); 316 } 317 } 318 } 319 320 PBM pbm(bpg); 321 pbm.run(); 322 checkMatching(bpg, pbm); 323 324 // do not always use the highest active node 325 PBM pbm2(bpg); 326 PBM::Elevator_t elev(bpg); 327 pbm2.elevator(elev); 328 pbm2.init(); 329 while (elev.highestActive() != INVALID) { 330 for (BlueNodeIt b(bpg); b!=INVALID; ++b) { 331 if (elev.active(b)) pbm2.push_relabel(b); 332 } 333 } 334 pbm2.writeMatching(); 335 checkMatching(bpg, pbm2); 336 } 337 338 339 return 0; 340 } 341 342 -
new file lemon/preflow_bp_matching_2.h
# HG changeset patch # User Daniel Poroszkai <poroszd@inf.elte.hu> # Date 1330129984 -3600 # Node ID ba89bdcf59f238ea9917f1d5d9515760c4847df7 # Parent b348df30553a63d6ecda56c1c2f1e1f18401f87a Other variant of preflow based bipartite matching diff --git a/lemon/preflow_bp_matching_2.h b/lemon/preflow_bp_matching_2.h new file mode 100644
- + 1 /* -*- mode: C++; indent-tabs-mode: nil; -*- 2 * 3 * This file is a part of LEMON, a generic C++ optimization library. 4 * 5 * Copyright (C) 2003-2012 6 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport 7 * (Egervary Research Group on Combinatorial Optimization, EGRES). 8 * 9 * Permission to use, modify and distribute this software is granted 10 * provided that this copyright notice appears in all copies. For 11 * precise terms see the accompanying LICENSE file. 12 * 13 * This software is provided "AS IS" with no warranty of any kind, 14 * express or implied, and with no claim as to its suitability for any 15 * purpose. 16 * 17 */ 18 19 #ifndef LEMON_PREFLOW_BP_MATCHING_2_H 20 #define LEMON_PREFLOW_BP_MATCHING_2_H 21 22 #include <lemon/core.h> 23 #include <lemon/elevator.h> 24 #include <vector> 25 #include <queue> 26 27 /// \ingroup matching 28 /// \file 29 /// \brief A preflow based bipartite matching algorithm. 30 namespace lemon { 31 32 /// \brief A preflow based bipartite matching algorithm 33 /// 34 /// Finds maximum cardinality matching in a given bipartite 35 /// graph using specialized preflow algorithm. 36 template <typename BPG> 37 class PreflowBpMatching2 { 38 public: 39 /// Type of the bipartite graph 40 typedef BPG BpGraph; 41 /// Type of the matching map 42 typedef typename BPG::template NodeMap<typename BPG::Edge> MatchingMap; 43 /// Type of the elevator. 44 typedef Elevator<BPG, typename BPG::BlueNode> Elevator_t; 45 46 private: 47 TEMPLATE_BPGRAPH_TYPEDEFS(BPG); 48 49 protected: 50 typedef typename BPG::template RedNodeMap<typename BPG::Arc> Preflow; 51 typedef typename BPG::template BlueNodeMap<int> Excess; 52 typedef std::vector<typename BPG::BlueNode> BlueNodes; 53 typedef std::vector<typename BPG::Arc> BlueFlow; 54 typedef typename BPG::template BlueNodeMap<BlueFlow> IncomingFlows; 55 56 const BpGraph& _bpg; 57 MatchingMap* _matching; 58 bool _local_matching; 59 Elevator_t* _elevator; 60 bool _local_elevator; 61 Preflow _flow; 62 Excess _excess; 63 64 void createStructures() { 65 if (!_matching) { 66 _matching = new MatchingMap(_bpg, INVALID); 67 _local_matching = true; 68 } 69 if (!_elevator) { 70 _elevator = new Elevator_t(_bpg, std::min(countRedNodes(_bpg) + 1, 71 countBlueNodes(_bpg))); 72 _local_elevator = true; 73 } 74 } 75 76 void destroyStructures() { 77 if (_local_matching) { 78 delete _matching; 79 } 80 if (_local_elevator) { 81 delete _elevator; 82 } 83 } 84 85 /// \brief Initialize the elevator 86 /// 87 /// When this function is called, the elevator is initialized with 88 /// the exact distance labels corresponding the current matching. 89 void initElevator() { 90 BlueNodes act_blue_nodes, next_blue_nodes; 91 BoolBlueNodeMap reached(_bpg, false); 92 IncomingFlows incoming(_bpg); 93 94 // Set the initial flow: 95 // choose arbitrary outgoing arc for every unmatched red node. 96 for (RedNodeIt r(_bpg); r!=INVALID; ++r) { 97 Arc a = (*_matching)[r] != INVALID ? 98 _bpg.direct((*_matching)[r], true) : OutArcIt(_bpg, r); 99 if (a != INVALID) { 100 _flow[r] = a; 101 BlueNode b = _bpg.asBlueNodeUnsafe(_bpg.target(a)); 102 incoming[b].push_back(a); 103 if (++_excess[b] == 2) { 104 act_blue_nodes.push_back(b); 105 reached.set(b, true); 106 } 107 } 108 } 109 // Now we start a Bfs from the blue nodes with excess 110 // (those which get flow from more than one red node) in the 111 // residual graph, to supply the exact distance labels to the elevator. 112 _elevator->initStart(); 113 while (!act_blue_nodes.empty()) { 114 for (typename BlueNodes::iterator b_it = act_blue_nodes.begin(); 115 b_it != act_blue_nodes.end(); ++b_it) { 116 117 _elevator->initAddItem(*b_it); 118 for (typename BlueFlow::iterator b_in = incoming[*b_it].begin(); 119 b_in != incoming[*b_it].end(); ++b_in) { 120 RedNode r = _bpg.asRedNodeUnsafe(_bpg.source(*b_in)); 121 for (OutArcIt r_oi(_bpg, r); r_oi!=INVALID; ++r_oi) { 122 BlueNode b = _bpg.asBlueNodeUnsafe(_bpg.target(r_oi)); 123 if (!reached[b]) { 124 next_blue_nodes.push_back(b); 125 reached.set(b, true); 126 } 127 } 128 } 129 } 130 _elevator->initNewLevel(); 131 next_blue_nodes.swap(act_blue_nodes); 132 next_blue_nodes.clear(); 133 } 134 _elevator->initFinish(); 135 136 // Activate blue nodes with 0 excess, i. e. those 137 // which does not get any flow. 138 for (BlueNodeIt b(_bpg); b != INVALID; ++b) { 139 if (_excess[b] == 0 && (*_elevator)[b] < _elevator->maxLevel()) { 140 _elevator->activate(b); 141 } 142 } 143 } 144 145 PreflowBpMatching2() {} 146 147 public: 148 149 /// \brief Constructor 150 /// 151 /// Constructs the class on the given bipartite graph. 152 PreflowBpMatching2(const BpGraph& bpg) : _bpg(bpg), 153 _matching(0), 154 _local_matching(false), 155 _elevator(0), 156 _local_elevator(false), 157 _flow(bpg, INVALID), 158 _excess(_bpg, 0) 159 160 {} 161 162 /// \brief Destructor 163 /// 164 /// Destructor. 165 ~PreflowBpMatching2() { 166 destroyStructures(); 167 } 168 169 /// \brief Sets the matching map 170 /// 171 /// Sets the matching map. It should contain a valid matching, 172 /// for example the empty matching is fine. 173 /// If you don't use this function before calling \ref init(), 174 /// an instance will be allocated automatically. 175 /// The destructor deallocates this automatically allocated map, 176 /// of course. 177 /// 178 /// \return <tt>(*this)</tt> 179 PreflowBpMatching2& matchingMap(MatchingMap& map) { 180 _matching = ↦ 181 _local_matching = false; 182 return *this; 183 } 184 185 /// \brief Returns a const reference to the matching map. 186 /// 187 /// Returns a const reference to the matching map, which contains 188 /// the matching edge for each node (or INVALID for unmatched nodes). 189 const MatchingMap& matchingMap() const { 190 return *_matching; 191 } 192 193 /// \brief Sets the elevator used by algorithm. 194 /// 195 /// Sets the elevator. 196 /// If you don't use this function before calling \ref init(), 197 /// an instance will be allocated automatically. 198 /// The destructor deallocates this automatically allocated elevator, 199 /// of course. 200 /// \return <tt>(*this)</tt> 201 PreflowBpMatching2& elevator(Elevator_t& elevator) { 202 _elevator = &elevator; 203 _local_elevator = false; 204 return *this; 205 } 206 207 /// \brief Returns a const reference to the elevator. 208 /// 209 /// Returns a const reference to the elevator, which keeps track 210 /// of the distance labels of nodes. 211 const Elevator_t& elevator() const { 212 return *_elevator; 213 } 214 215 /// \brief Initializes the algorithm 216 /// 217 /// Allocates the elevator and the matching map if necessary, 218 /// and initializes the elevator. 219 /// 220 /// \pre \ref init() is not called. 221 void init() { 222 createStructures(); 223 initElevator(); 224 } 225 226 /// \brief Smarter initialization of the matching. 227 /// 228 /// Allocates the matching map if necessary, and initializes 229 /// the algorithm with a trivially found matching: iterating 230 /// on every edge, take it in if both endnodes are unmatched. 231 /// 232 /// \return Size of the so found matching. 233 /// 234 /// \note heuristicInit() replaces init() regarding the preconditions. 235 int heuristicInit() { 236 createStructures(); 237 int size = 0; 238 239 for (BlueNodeIt b(_bpg); b!=INVALID; ++b) { 240 bool matched = false; 241 for (OutArcIt oi(_bpg, b); oi != INVALID && !matched; ++oi) { 242 if ((*_matching)[_bpg.target(oi)] == INVALID) { 243 _matching->set(_bpg.source(oi), oi); 244 _matching->set(_bpg.target(oi), oi); 245 matched = true; 246 ++size; 247 } 248 } 249 } 250 251 initElevator(); 252 return size; 253 } 254 255 /// \brief Initialize the matching from a map. 256 /// 257 /// Allocates the matching map if necessary, and 258 /// initializes the algorithm with a matching given as a bool edgemap, 259 /// in which an edge is true if it is in the matching. 260 /// Also initializes the elevator. 261 /// 262 /// \return \c true if the matching is valid. 263 /// 264 /// \note matchingInit() replaces init() regarding the preconditions. 265 bool matchingInit(const BoolEdgeMap& matching) { 266 createStructures(); 267 268 bool valid = true; 269 for(EdgeIt it(_bpg); it!=INVALID; ++it) { 270 if (matching[it]) { 271 Node red = _bpg.redNode(it); 272 Node blue = _bpg.blueNode(it); 273 if ((*_matching)[red] != INVALID || (*_matching)[blue] != INVALID) { 274 valid = false; 275 } else { 276 _matching->set(red, it); 277 _matching->set(blue, it); 278 } 279 } 280 } 281 282 initElevator(); 283 return valid; 284 } 285 286 287 /// \brief Perform a pull/relabel operation on a node 288 /// 289 /// Perform a pull/relabel operation on a node: 290 /// search an incident blue node on one level lower, 291 /// and flip its flow to the 'chosen', or relabel 292 /// if none found. 293 /// 294 /// \pre \ref init() is called, and the node is active. 295 void pull_relabel(const BlueNode& chosen) { 296 297 int min_inc_level = _elevator->maxLevel() + 1; 298 for (OutArcIt b_oi(_bpg, chosen); b_oi!=INVALID; ++b_oi) { 299 RedNode r = _bpg.asRedNodeUnsafe(_bpg.target(b_oi)); 300 BlueNode lower = _bpg.asBlueNodeUnsafe(_bpg.target(_flow[r])); 301 if ((*_elevator)[(lower)] < (*_elevator)[chosen]) { 302 if (--_excess[lower] == 0) { 303 _elevator->activate(lower); 304 } 305 _flow[r] = _bpg.direct(b_oi, true); 306 ++_excess[chosen]; 307 _elevator->deactivate(chosen); 308 break; 309 } else { 310 if ((*_elevator)[lower] < min_inc_level && 311 lower != chosen) { 312 min_inc_level = (*_elevator)[lower]; 313 } 314 } 315 } 316 // If none found, then relabel. 317 if (_excess[chosen] == 0) { 318 int chosen_level = (*_elevator)[chosen]; 319 if (min_inc_level >= _elevator->maxLevel() - 1) { 320 _elevator->lift(chosen, _elevator->maxLevel()); 321 _elevator->deactivate(chosen); 322 } else { 323 _elevator->lift(chosen, min_inc_level+1); 324 } 325 if (_elevator->emptyLevel(chosen_level)) { 326 _elevator->liftToTop(chosen_level); 327 } 328 } 329 } 330 331 /// \brief Executes the algorithm 332 /// 333 /// Executes pull/relabel operation on the highest active node 334 /// until there is. After start(), writeMatching() should be called to 335 /// remove unnecessary flow and choose a valid maximum matching. 336 /// 337 /// \pre \ref init() must be called before using this function. 338 void start() { 339 while (_elevator->highestActive() != INVALID) 340 pull_relabel(_elevator->highestActive()); 341 } 342 343 /// \brief Deduces and writes a matching from the preflow 344 /// 345 /// The underlying preflow algorithm does not use the matching map, 346 /// so before using the matching query functions ( \ref mate(), \ref 347 /// matchingSize(), etc.), you should explicitly call this function 348 /// to write a matching based on the current preflow to the matching 349 /// map. When the preflow results an ambiguous matching (a node has 350 /// incoming flow on more than one arc), it chooses arbitrary one. 351 void writeMatching() { 352 for (NodeIt it(_bpg); it != INVALID; ++it) { 353 _matching->set(it, INVALID); 354 } 355 for (RedNodeIt r(_bpg); r != INVALID; ++r) { 356 if (_flow[r] != INVALID) { 357 BlueNode mate = _bpg.blueNode(_flow[r]); 358 if (_excess[mate] > 1) { 359 --_excess[mate]; 360 } else { 361 _matching->set(mate, _flow[r]); 362 _matching->set(r, _flow[r]); 363 } 364 } 365 } 366 } 367 368 /// \brief Runs the algorithm. 369 /// 370 /// bp1.run() is just a shorthand for: 371 /// 372 ///\code 373 /// bp1.heuristicInit(); 374 /// bp1.start(); 375 /// bp1.writeMatching(); 376 ///\endcode 377 void run() { 378 heuristicInit(); 379 start(); 380 writeMatching(); 381 } 382 383 /// \brief Size of the matching 384 /// 385 /// Returns the size of the current matching. 386 /// Function runs in \f$O(n)\f$ time. 387 int matchingSize() const { 388 int size = 0; 389 for (RedNodeIt it(_bpg); it!=INVALID; ++it) { 390 if ((*_matching)[it] != INVALID) ++size; 391 } 392 return size; 393 } 394 395 /// \brief Return \c true if the given edge is in the matching. 396 /// 397 /// This function returns \c true if the given edge is in the current 398 /// matching. 399 bool matching(const Edge& edge) const { 400 return edge == (*_matching)[_bpg.redNode(edge)]; 401 } 402 403 /// \brief Return the matching edge incident to the given node. 404 /// 405 /// This function returns the matching edge incident to the 406 /// given node in the current matching or \c INVALID if the node is 407 /// not covered by the matching. 408 Edge matching(const Node& n) const { 409 return (*_matching)[n]; 410 } 411 412 /// \brief Return the mate of the given node. 413 /// 414 /// This function returns the mate of the given node in the current 415 /// matching or \c INVALID if the node is not covered by the matching. 416 Node mate(const Node& n) const { 417 return (*_matching)[n] != INVALID ? 418 _bpg.oppositeNode(n, (*_matching)[n]) : INVALID; 419 } 420 421 /// \brief Query a vertex cover. 422 /// 423 /// This function finds a vertex cover based on the current matching, 424 /// and in \c cover sets node in the cover \c true and \c false 425 /// otherwise. In case of maximal matching, this will be a minimal 426 /// vertex cover. 427 /// Function runs in \f$O(n + e)\f$ time. 428 /// 429 /// \return The number of nodes in the cover. 430 int cover(BoolNodeMap &cover) const { 431 int count = 0; 432 std::queue<RedNode> q; 433 434 for (BlueNodeIt b(_bpg); b!=INVALID; ++b) 435 cover[b] = false; 436 437 for (RedNodeIt r(_bpg); r!=INVALID; ++r) { 438 if ((*_matching)[r] == INVALID) { 439 cover[r] = false; 440 q.push(r); 441 } else { 442 cover[r] = true; 443 ++count; 444 } 445 } 446 447 while (!q.empty()) { 448 for (OutArcIt it(_bpg, q.front()); it!=INVALID; ++it) { 449 BlueNode blue(_bpg.asBlueNodeUnsafe(_bpg.target(it))); 450 if (cover[blue]) continue; 451 452 cover[blue] = true; 453 ++count; 454 455 if ((*_matching)[blue] != INVALID) { 456 RedNode nr = _bpg.redNode((*_matching)[blue]); 457 q.push(nr); 458 cover[nr] = false; 459 --count; 460 } 461 } 462 q.pop(); 463 } 464 465 return count; 466 } 467 468 /// \brief Query a barrier of red nodes. 469 /// 470 /// This function tries to create a barrier of red nodes, which should: 471 /// - contain every unmatched red nodes, 472 /// - contain exactly that many matched red nodes as much blue nodes 473 /// there are in its neighbor set. 474 /// 475 /// Barrier exist if and only if the matching is maximal. If exist, the 476 /// map is set true for nodes of the barrier, and false for the others. 477 /// 478 /// Function runs in \f$O(n + e)\f$ time. 479 /// 480 /// \return \c true if a barrier found. 481 bool redBarrier(BoolRedNodeMap &barrier) const { 482 bool exist = true; 483 std::queue<RedNode> q; 484 485 for (RedNodeIt r(_bpg); r!=INVALID; ++r) { 486 if ((*_matching)[r] == INVALID) { 487 barrier[r] = true; 488 q.push(r); 489 } else { 490 barrier[r] = false; 491 } 492 } 493 494 while (!q.empty() && exist) { 495 for (OutArcIt it(_bpg, q.front()); it!=INVALID; ++it) { 496 BlueNode blue(_bpg.asBlueNodeUnsafe(_bpg.target(it))); 497 498 if (exist &= ((*_matching)[blue] != INVALID)) { 499 RedNode r = _bpg.redNode((*_matching)[blue]); 500 if (!barrier[r]) { 501 q.push(r); 502 barrier[r] = true; 503 } 504 } 505 } 506 q.pop(); 507 } 508 509 return exist; 510 } 511 512 /// \brief Query a barrier of blue nodes. 513 /// 514 /// This function tries to create a barrier of blue nodes, which should: 515 /// - contain every unmatched blue nodes, 516 /// - contain exactly that many matched blue nodes as much red nodes 517 /// there are in its neighbor set. 518 /// 519 /// Barrier exist if and only if the matching is maximal. If exist, the 520 /// map is set true for nodes of the barrier, and false for the others. 521 /// 522 /// Function runs in \f$O(n + e)\f$ time. 523 /// 524 /// \return \c true if a barrier found. 525 bool blueBarrier(BoolBlueNodeMap &barrier) const { 526 bool exist = true; 527 std::queue<BlueNode> q; 528 529 for (BlueNodeIt b(_bpg); b!=INVALID; ++b) { 530 if ((*_matching)[b] == INVALID) { 531 barrier[b] = true; 532 q.push(b); 533 } else { 534 barrier[b] = false; 535 } 536 } 537 538 while (!q.empty() && exist) { 539 for (OutArcIt it(_bpg, q.front()); it!=INVALID; ++it) { 540 RedNode r(_bpg.asRedNodeUnsafe(_bpg.target(it))); 541 542 if (exist &= ((*_matching)[r] != INVALID)) { 543 BlueNode b = _bpg.blueNode((*_matching)[r]); 544 if (!barrier[b]) { 545 q.push(b); 546 barrier[b] = true; 547 } 548 } 549 } 550 q.pop(); 551 } 552 553 return exist; 554 } 555 }; 556 557 } 558 #endif -
test/CMakeLists.txt
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
a b 48 48 path_test 49 49 planarity_test 50 50 preflow_test 51 preflow_bp_matching_2_test 51 52 preflow_bp_matching_test 52 53 radix_sort_test 53 54 random_test -
new file test/preflow_bp_matching_2_test.cc
diff --git a/test/preflow_bp_matching_2_test.cc b/test/preflow_bp_matching_2_test.cc new file mode 100644
- + 1 /* -*- mode: C++; indent-tabs-mode: nil; -*- 2 * 3 * This file is a part of LEMON, a generic C++ optimization library. 4 * 5 * Copyright (C) 2003-2012 6 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport 7 * (Egervary Research Group on Combinatorial Optimization, EGRES). 8 * 9 * Permission to use, modify and distribute this software is granted 10 * provided that this copyright notice appears in all copies. For 11 * precise terms see the accompanying LICENSE file. 12 * 13 * This software is provided "AS IS" with no warranty of any kind, 14 * express or implied, and with no claim as to its suitability for any 15 * purpose. 16 * 17 */ 18 19 #include <iostream> 20 #include <sstream> 21 #include <list> 22 23 #include <lemon/random.h> 24 #include <lemon/preflow_bp_matching_2.h> 25 #include <lemon/smart_graph.h> 26 #include <lemon/concepts/bpgraph.h> 27 #include <lemon/concepts/maps.h> 28 #include <lemon/lgf_reader.h> 29 30 #include "test_tools.h" 31 32 using namespace std; 33 using namespace lemon; 34 35 const int lgfn = 3; 36 const string lgf[lgfn] = { 37 // Empty graph 38 "@red_nodes\n" 39 "label\n" 40 "\n" 41 "@blue_nodes\n" 42 "label\n" 43 "\n" 44 "@edges\n" 45 " label\n" 46 "\n", 47 48 // No red nodes 49 "@red_nodes\n" 50 "label\n" 51 "\n" 52 "@blue_nodes\n" 53 "label\n" 54 "1\n" 55 "@edges\n" 56 " label\n" 57 "\n", 58 59 // No blue nodes 60 "@red_nodes\n" 61 "label\n" 62 "1\n" 63 "@blue_nodes\n" 64 "label\n" 65 "\n" 66 "@edges\n" 67 " label\n" 68 "\n", 69 }; 70 71 const string u_lgf = 72 "@red_nodes\n" 73 "label\n" 74 "1\n" 75 "2\n" 76 "3\n" 77 "@blue_nodes\n" 78 "label\n" 79 "4\n" 80 "5\n" 81 "6\n" 82 "7\n" 83 "@edges\n" 84 "label\n" 85 "1 4 1\n" 86 "1 5 2\n" 87 "2 5 3\n" 88 "2 6 4\n" 89 "2 7 5\n" 90 "3 7 6\n" 91 ; 92 93 94 95 void checkPreflowBpMatching2Compile() { 96 typedef concepts::BpGraph BpGraph; 97 BPGRAPH_TYPEDEFS(BpGraph); 98 typedef PreflowBpMatching2<BpGraph> PBM; 99 100 BpGraph graph; 101 RedNode red; 102 BlueNode blue; 103 Node node; 104 Edge edge; 105 BoolEdgeMap init_matching(graph); 106 PBM pbm_test(graph); 107 const PBM& const_pbm_test = pbm_test; 108 PBM::MatchingMap matching_map(graph, INVALID); 109 PBM::Elevator_t elevator(graph); 110 BpGraph::NodeMap<bool> cov(graph); 111 int cover_size; 112 BoolRedNodeMap rb(graph); 113 BoolBlueNodeMap bb(graph); 114 115 pbm_test.matchingMap(matching_map); 116 pbm_test.elevator(elevator); 117 pbm_test.init(); 118 pbm_test.heuristicInit(); 119 pbm_test.matchingInit(init_matching); 120 121 pbm_test.start(); 122 pbm_test.pull_relabel(blue); 123 pbm_test.run(); 124 pbm_test.writeMatching(); 125 126 const_pbm_test.matchingSize(); 127 const_pbm_test.matching(node); 128 const_pbm_test.matching(edge); 129 const_pbm_test.mate(red); 130 const_pbm_test.mate(blue); 131 const_pbm_test.mate(node); 132 const PBM::MatchingMap& max_matching = const_pbm_test.matchingMap(); 133 edge = max_matching[node]; 134 const PBM::Elevator_t& elev = const_pbm_test.elevator(); 135 blue = elev.highestActive(); 136 cover_size = const_pbm_test.cover(cov); 137 const_pbm_test.redBarrier(rb); 138 const_pbm_test.blueBarrier(bb); 139 } 140 141 typedef SmartBpGraph BpGraph; 142 typedef PreflowBpMatching2<BpGraph> PBM; 143 BPGRAPH_TYPEDEFS(BpGraph); 144 145 void checkMatchingValidity(const BpGraph& bpg, const PBM& pbm) { 146 BoolBlueNodeMap matched(bpg, false); 147 for (RedNodeIt it(bpg); it != INVALID; ++it) { 148 if (pbm.matching(it) != INVALID) { 149 check(pbm.mate(pbm.mate(it)) == it, "Wrong matching"); 150 matched.set(bpg.asBlueNode(pbm.mate(it)), true); 151 } 152 } 153 for (BlueNodeIt it(bpg); it != INVALID; ++it) { 154 // != is used as logical xor here 155 check(matched[it] != (pbm.matching(it) == INVALID), "Wrong matching"); 156 } 157 } 158 159 void checkCover(const BpGraph& bpg, const PBM& pbm) { 160 int cov_size; 161 BoolNodeMap cov(bpg); 162 cov_size = pbm.cover(cov); 163 BoolEdgeMap covered(bpg, false); 164 for (NodeIt n(bpg); n!=INVALID; ++n) { 165 for (IncEdgeIt e(bpg, n); e!=INVALID; ++e) { 166 covered[e] = true; 167 } 168 } 169 bool all = true; 170 for (EdgeIt e(bpg); e!=INVALID; ++e) { 171 all&=covered[e]; 172 } 173 174 check(all, "Invalid cover"); 175 check(pbm.matchingSize() == cov_size, "Suboptimal matching."); 176 } 177 void checkBarriers(const BpGraph& bpg, const PBM& pbm) { 178 BoolRedNodeMap rb(bpg); 179 BoolBlueNodeMap bb(bpg); 180 181 check(pbm.redBarrier(rb), "Suboptimal matching."); 182 check(pbm.blueBarrier(bb), "Suboptimal matching."); 183 184 BoolNodeMap reached(bpg, false); 185 186 for (RedNodeIt r(bpg); r!=INVALID; ++r) { 187 if (!rb[r]) continue; 188 for (OutArcIt a(bpg, r); a!=INVALID; ++a) { 189 reached.set(bpg.target(a), true); 190 } 191 } 192 for (BlueNodeIt b(bpg); b!=INVALID; ++b) { 193 if (!bb[b]) continue; 194 for (OutArcIt a(bpg, b); a!=INVALID; ++a) { 195 reached.set(bpg.target(a), true); 196 } 197 } 198 199 int ur = 0, ub = 0, rc = 0, bc = 0; 200 for (RedNodeIt r(bpg); r!=INVALID; ++r) { 201 if (rb[r]) ++rc; 202 if (pbm.mate(r) == INVALID) ++ur; 203 if (reached[r]) --bc; 204 } 205 206 for (BlueNodeIt b(bpg); b!=INVALID; ++b) { 207 if (bb[b]) ++bc; 208 if (pbm.mate(b) == INVALID) ++ub; 209 if (reached[b]) --rc; 210 } 211 212 check(rc == ur, "Red barrier is wrong."); 213 check(bc == ub, "Blue barrier is wrong."); 214 } 215 216 void checkMatching(const BpGraph& bpg, const PBM& pbm) { 217 checkMatchingValidity(bpg, pbm); 218 checkBarriers(bpg, pbm); 219 checkCover(bpg, pbm); 220 } 221 222 223 int main() { 224 // Checking hard wired extremal graphs 225 for (int i=0; i<lgfn; ++i) { 226 BpGraph bpg; 227 istringstream lgfs(lgf[i]); 228 bpGraphReader(bpg, lgfs).run(); 229 PBM pbm(bpg); 230 pbm.run(); 231 checkMatching(bpg, pbm); 232 } 233 234 // Checking some use cases 235 BpGraph bpg; 236 istringstream u_lgfs(u_lgf); 237 bpGraphReader(bpg, u_lgfs).run(); 238 239 { 240 PBM pbm(bpg); 241 pbm.init(); 242 BpGraph::BlueNode b = pbm.elevator().highestActive(); 243 if (b != INVALID) { 244 pbm.pull_relabel(b); 245 } 246 pbm.start(); 247 pbm.writeMatching(); 248 check(pbm.elevator().highestActive() == INVALID, "Error in elevator"); 249 checkMatching(bpg, pbm); 250 } 251 252 { 253 PBM pbm(bpg); 254 pbm.heuristicInit(); 255 BpGraph::BlueNode b = pbm.elevator().highestActive(); 256 if (b != INVALID) { 257 pbm.pull_relabel(b); 258 } 259 pbm.start(); 260 checkMatching(bpg, pbm); 261 } 262 { 263 PBM pbm(bpg); 264 PBM::MatchingMap m(bpg, INVALID); 265 pbm.matchingMap(m); 266 pbm.run(); 267 checkMatching(bpg, pbm); 268 269 } 270 { 271 PBM pbm(bpg); 272 BpGraph::EdgeMap<bool> init_matching(bpg, false); 273 BpGraph::Edge e; 274 bpg.first(e); 275 init_matching[e] = true; 276 check(pbm.matchingInit(init_matching), 277 "Wrong result from matchingInit()"); 278 pbm.start(); 279 pbm.writeMatching(); 280 checkMatching(bpg, pbm); 281 } 282 { 283 PBM pbm(bpg); 284 PBM::MatchingMap m(bpg, INVALID); 285 pbm.matchingMap(m); 286 BpGraph::EdgeMap<bool> init_matching(bpg, true); 287 check(!pbm.matchingInit(init_matching), 288 "Wrong result from matchingInit()"); 289 pbm.start(); 290 pbm.writeMatching(); 291 checkMatching(bpg, pbm); 292 } 293 294 // Checking some random generated graphs 295 const int random_test_n = 10; 296 const int max_red_n = 1000; 297 const int min_red_n = 10; 298 const int max_blue_n = 1000; 299 const int min_blue_n = 10; 300 for (int i=0; i<random_test_n; ++i) { 301 int red_n = rnd.integer(min_red_n, max_red_n); 302 int blue_n = rnd.integer(min_blue_n, max_blue_n); 303 double prob = rnd(); 304 305 SmartBpGraph bpg; 306 for (int j=0; j<red_n; ++j) { 307 bpg.addRedNode(); 308 } 309 for (int j=0; j<blue_n; ++j) { 310 bpg.addBlueNode(); 311 } 312 for (SmartBpGraph::RedNodeIt r(bpg); r!=INVALID; ++r) { 313 for (SmartBpGraph::BlueNodeIt b(bpg); b!=INVALID; ++b) { 314 if (rnd() < prob/10.) { 315 bpg.addEdge(r,b); 316 } 317 } 318 } 319 320 PBM pbm(bpg); 321 pbm.run(); 322 checkMatching(bpg, pbm); 323 324 // do not always use the highest active node 325 PBM pbm2(bpg); 326 PBM::Elevator_t elev(bpg); 327 pbm2.elevator(elev); 328 pbm2.init(); 329 while (elev.highestActive() != INVALID) { 330 for (BlueNodeIt b(bpg); b!=INVALID; ++b) { 331 if (elev.active(b)) pbm2.pull_relabel(b); 332 } 333 } 334 pbm2.writeMatching(); 335 checkMatching(bpg, pbm2); 336 } 337 338 339 return 0; 340 } 341 342 -
tools/bp-matching-benchmark.cc
# HG changeset patch # User Daniel Poroszkai <poroszd@inf.elte.hu> # Date 1329867694 -3600 # Node ID 93fa3bb078dd0fb7e198d73c42999ceccd1bbe27 # Parent add2bdd7f7a635e1f7d080dd433d342d9271e779 Benchmark tool extended diff --git a/tools/bp-matching-benchmark.cc b/tools/bp-matching-benchmark.cc
a b 23 23 #include <lemon/hopcroft_karp.h> 24 24 #include <lemon/matching.h> 25 25 #include <lemon/preflow.h> 26 #include <lemon/preflow_bp_matching.h> 27 #include <lemon/preflow_bp_matching_2.h> 26 28 #include <lemon/bp_matching.h> 27 29 #include <lemon/cost_scaling.h> 28 30 #include <lemon/capacity_scaling.h> … … 32 34 using namespace std; 33 35 using namespace lemon; 34 36 37 void check(bool c, const string &msg) { 38 if (!c) { 39 std::cerr << "Error: " << msg << "\n"; 40 } 41 } 35 42 36 43 template <typename BpGraph, typename Digraph> 37 44 void create_network(const BpGraph &bpg, … … 108 115 /// 109 116 /// Depending on the input, it runs the \ref lemon::HopcroftKarp 110 117 /// "Hopcroft-Karp", the \ref lemon::MaxMatching "general matching" 111 /// and the \ref lemon::Preflow "Preflow" algorithms for finding 118 /// and three algorithms based on preflow (\ref lemon::Preflow "preflow 119 /// on extended graph", \ref lemon::PreflowBpMatching, 120 /// \ref lemon::PreflowBpMatching2") for finding 112 121 /// a maximum cardinality bipartite matching, or the 113 122 /// \ref lemon::MaxWeightedBpMatching "max. weighted matching 114 123 /// for sparse bipartite graphs", the \ref lemon::MaxWeightedMatching … … 118 127 /// simplex" algorithms to find the maximum weighted matching. 119 128 int main() { 120 129 DimacsDescriptor desc = dimacsType(cin); 121 int red_n, blue_n, edge_n = desc.edgeNum;122 130 if (desc.type == DimacsDescriptor::UBM) { 123 131 // unweighted bipartite matching 124 Timer hk_t(false), mm_t(false), pf_t(false); 125 int hk_r, mm_r, pf_r; 132 Timer hk_t(false), 133 mm_t(false), 134 pf_t(false), 135 pbm_t(false), 136 pbm2_t(false); 126 137 127 138 BpGraph *bpg = new BpGraph(); 128 139 readDimacsUbm(cin, *bpg, desc); 129 red_n = countRedNodes(*bpg);130 blue_n = countBlueNodes(*bpg);131 140 132 141 { 133 142 hk_t.start(); 134 143 HopcroftKarp<BpGraph> hk(*bpg); 135 144 hk.run(); 136 hk_r = hk.matchingSize();137 145 hk_t.stop(); 146 BpGraph::NodeMap<bool> m(*bpg); 147 check(hk.matchingSize() == hk.cover(m), 148 "Suboptimal matching in Hopcroft-Karp."); 138 149 } 139 150 { 140 151 mm_t.start(); 141 152 MaxMatching<BpGraph> mm(*bpg); 142 153 mm.run(); 143 mm_r = mm.matchingSize();144 154 mm_t.stop(); 145 155 } 156 { 157 pbm_t.start(); 158 PreflowBpMatching<BpGraph> pbm(*bpg); 159 pbm.run(); 160 pbm_t.stop(); 161 BpGraph::NodeMap<bool> m(*bpg); 162 check(pbm.matchingSize() == pbm.cover(m), 163 "Suboptimal matching in push-relabel."); 164 } 165 { 166 pbm2_t.start(); 167 PreflowBpMatching2<BpGraph> pbm2(*bpg); 168 pbm2.run(); 169 pbm2_t.stop(); 170 BpGraph::NodeMap<bool> m(*bpg); 171 check(pbm2.matchingSize() == pbm2.cover(m), 172 "Suboptimal matching in pull-relabel."); 173 } 146 174 147 175 Digraph *network = new Digraph(); 148 176 Digraph::Node source, target; 149 177 create_network(*bpg, *network, source, target); 150 delete bpg;151 178 { 152 179 pf_t.start(); 153 180 Preflow<Digraph, ConstMap<Digraph::Arc, long> > 154 181 pf(*network, constMap<Digraph::Arc, long>(1), source, target); 155 182 pf.run(); 156 pf_r = pf.flowValue();157 183 pf_t.stop(); 158 184 } 159 185 delete network; 160 186 161 187 cout << "Benchmarking in unweighted case, using a bipartite graph\n" 162 << "with " << red_n << " red, " << blue_n << " blue nodes, and " 163 << edge_n << " edges.\n" 188 << "with " << countRedNodes(*bpg) << " red, " 189 << countBlueNodes(*bpg) << " blue nodes, and " 190 << desc.edgeNum << " edges.\n" 164 191 << "--------------------------------------------------------\n" 165 << "Algorithm used Matching size Time\n" 166 << "Hopcroft-Karp " << hk_r 167 << " " << hk_t.realTime() << "\n" 168 << "General matching " << mm_r 169 << " " << mm_t.realTime() << "\n" 170 << "Preflow " << pf_r 171 << " " << pf_t.realTime() << endl; 192 << "Algorithm used Time\n" 193 << "Hopcroft-Karp " << hk_t.realTime() << "\n" 194 << "General matching " << mm_t.realTime() << "\n" 195 << "Push-relabel matching " << pbm_t.realTime() << "\n" 196 << "Pull-relabel matching " << pbm2_t.realTime() << "\n" 197 << "Preflow " << pf_t.realTime() << endl; 172 198 199 delete bpg; 173 200 } else if (desc.type == DimacsDescriptor::WBM) { 174 201 // weighted bipartite matching 175 202 Timer bm_t(false), … … 182 209 BpGraph *bpg = new BpGraph(); 183 210 BpGraph::EdgeMap<long> *weight = new BpGraph::EdgeMap<long>(*bpg); 184 211 readDimacsWbm(cin, *bpg, *weight, desc); 185 red_n = countRedNodes(*bpg);186 blue_n = countBlueNodes(*bpg);187 212 188 213 { 189 214 bm_t.start(); … … 241 266 delete supply; 242 267 243 268 cout << "Benchmarking on weighted bipartite matching problem on a " 244 << "graph\nwith " << red_n << " red, " << blue_n 245 << " blue nodes and " << edge_n << " edges\n" 269 << "graph\nwith " << countRedNodes(*bpg) << " red, " 270 << countBlueNodes(*bpg) << " blue nodes and " 271 << desc.edgeNum << " edges\n" 246 272 << "--------------------------------------------------------\n" 247 273 << "Algorithm used Maximum weight Time\n" 248 274 << "Bipartite matching " << bm_r