rustworkx_core/max_weight_matching.rs
1// Licensed under the Apache License, Version 2.0 (the "License"); you may
2// not use this file except in compliance with the License. You may obtain
3// a copy of the License at
4//
5// http://www.apache.org/licenses/LICENSE-2.0
6//
7// Unless required by applicable law or agreed to in writing, software
8// distributed under the License is distributed on an "AS IS" BASIS, WITHOUT
9// WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the
10// License for the specific language governing permissions and limitations
11// under the License.
12
13// Needed to pass shared state between functions
14// closures don't work because of recursion
15#![allow(clippy::too_many_arguments)]
16// Allow single character names to match naming convention from
17// paper
18#![allow(clippy::many_single_char_names)]
19
20use crate::dictmap::*;
21use std::cmp::max;
22use std::hash::Hash;
23use std::mem;
24
25use hashbrown::{HashMap, HashSet};
26
27use petgraph::Undirected;
28use petgraph::visit::{
29 EdgeCount, EdgeRef, GraphBase, GraphProp, IntoEdges, IntoNodeIdentifiers, NodeCount,
30 NodeIndexable,
31};
32
33/// Return 2 * slack of edge k (does not work inside blossoms).
34fn slack(edge_index: usize, dual_var: &[i128], edges: &[(usize, usize, i128)]) -> i128 {
35 let (source_index, target_index, weight) = edges[edge_index];
36 dual_var[source_index] + dual_var[target_index] - 2 * weight
37}
38
39/// Generate the leaf vertices of a blossom.
40fn blossom_leaves<E>(
41 blossom: usize,
42 num_nodes: usize,
43 blossom_children: &[Vec<usize>],
44) -> Result<Vec<usize>, E> {
45 let mut out_vec: Vec<usize> = Vec::new();
46 if blossom < num_nodes {
47 out_vec.push(blossom);
48 } else {
49 let child_blossom = &blossom_children[blossom];
50 for c in child_blossom {
51 let child = *c;
52 if child < num_nodes {
53 out_vec.push(child);
54 } else {
55 for v in blossom_leaves(child, num_nodes, blossom_children)? {
56 out_vec.push(v);
57 }
58 }
59 }
60 }
61 Ok(out_vec)
62}
63
64/// Assign label t to the top-level blossom containing vertex w
65/// and record the fact that w was reached through the edge with
66/// remote endpoint p.
67fn assign_label<E>(
68 w: usize,
69 t: usize,
70 p: Option<usize>,
71 num_nodes: usize,
72 in_blossoms: &[usize],
73 labels: &mut Vec<Option<usize>>,
74 label_ends: &mut Vec<Option<usize>>,
75 best_edge: &mut Vec<Option<usize>>,
76 queue: &mut Vec<usize>,
77 blossom_children: &[Vec<usize>],
78 blossom_base: &[Option<usize>],
79 endpoints: &[usize],
80 mate: &HashMap<usize, usize>,
81) -> Result<(), E> {
82 let b = in_blossoms[w];
83 assert!(labels[w] == Some(0) && labels[b] == Some(0));
84 labels[w] = Some(t);
85 labels[b] = Some(t);
86 label_ends[b] = p;
87 label_ends[w] = p;
88 best_edge[w] = None;
89 best_edge[b] = None;
90 if t == 1 {
91 // b became an S-vertex/blossom; add it(s vertices) to the queue
92 queue.append(&mut blossom_leaves(b, num_nodes, blossom_children)?);
93 } else if t == 2 {
94 // b became a T-vertex/blossom; assign label S to its mate.
95 // (If b is a non-trivial blossom, its base is the only vertex
96 // with an external mate.)
97 let blossom_index: usize = b;
98 let base: usize = blossom_base[blossom_index].unwrap();
99 assert!(mate.get(&base).is_some());
100 assign_label(
101 endpoints[mate[&base]],
102 1,
103 mate.get(&base).map(|p| p ^ 1),
104 num_nodes,
105 in_blossoms,
106 labels,
107 label_ends,
108 best_edge,
109 queue,
110 blossom_children,
111 blossom_base,
112 endpoints,
113 mate,
114 )?;
115 }
116 Ok(())
117}
118
119/// Trace back from vertices v and w to discover either a new blossom
120/// or an augmenting path. Return the base vertex of the new blossom or None.
121fn scan_blossom(
122 node_a: usize,
123 node_b: usize,
124 in_blossoms: &[usize],
125 blossom_base: &[Option<usize>],
126 endpoints: &[usize],
127 labels: &mut [Option<usize>],
128 label_ends: &[Option<usize>],
129 mate: &HashMap<usize, usize>,
130) -> Option<usize> {
131 let mut v: Option<usize> = Some(node_a);
132 let mut w: Option<usize> = Some(node_b);
133 // Trace back from v and w, placing breadcrumbs as we go
134 let mut path: Vec<usize> = Vec::new();
135 let mut base: Option<usize> = None;
136 while v.is_some() || w.is_some() {
137 // Look for a breadcrumb in v's blossom or put a new breadcrump
138 let mut blossom = in_blossoms[v.unwrap()];
139 if labels[blossom].is_none() || labels[blossom].unwrap() & 4 > 0 {
140 base = blossom_base[blossom];
141 break;
142 }
143 assert!(labels[blossom] == Some(1));
144 path.push(blossom);
145 labels[blossom] = Some(5);
146 // Trace one step back.
147 assert!(label_ends[blossom] == mate.get(&blossom_base[blossom].unwrap()).copied());
148 if label_ends[blossom].is_none() {
149 // The base of blossom is single; stop tracing this path
150 v = None;
151 } else {
152 let tmp = endpoints[label_ends[blossom].unwrap()];
153 blossom = in_blossoms[tmp];
154 assert!(labels[blossom] == Some(2));
155 // blossom is a T-blossom; trace one more step back.
156 assert!(label_ends[blossom].is_some());
157 v = Some(endpoints[label_ends[blossom].unwrap()]);
158 }
159 // Swap v and w so that we alternate between both paths.
160 if w.is_some() {
161 mem::swap(&mut v, &mut w);
162 }
163 }
164 // Remove breadcrumbs.
165 for blossom in path {
166 labels[blossom] = Some(1);
167 }
168 // Return base vertex, if we found one.
169 base
170}
171
172/// Construct a new blossom with given base, containing edge k which
173/// connects a pair of S vertices. Label the new blossom as S; set its dual
174/// variable to zero; relabel its T-vertices to S and add them to the queue.
175fn add_blossom<E>(
176 base: usize,
177 edge: usize,
178 blossom_children: &mut [Vec<usize>],
179 num_nodes: usize,
180 edges: &[(usize, usize, i128)],
181 in_blossoms: &mut [usize],
182 dual_var: &mut [i128],
183 labels: &mut [Option<usize>],
184 label_ends: &mut [Option<usize>],
185 best_edge: &mut [Option<usize>],
186 queue: &mut Vec<usize>,
187 blossom_base: &mut [Option<usize>],
188 endpoints: &[usize],
189 blossom_endpoints: &mut [Vec<usize>],
190 unused_blossoms: &mut Vec<usize>,
191 blossom_best_edges: &mut [Vec<usize>],
192 blossom_parents: &mut [Option<usize>],
193 neighbor_endpoints: &[Vec<usize>],
194 mate: &HashMap<usize, usize>,
195) -> Result<(), E> {
196 let (mut v, mut w, _weight) = edges[edge];
197 let blossom_b = in_blossoms[base];
198 let mut blossom_v = in_blossoms[v];
199 let mut blossom_w = in_blossoms[w];
200 // Create blossom
201 let blossom = unused_blossoms.pop().unwrap();
202 blossom_base[blossom] = Some(base);
203 blossom_parents[blossom] = None;
204 blossom_parents[blossom_b] = Some(blossom);
205 // Make list of sub-blossoms and their interconnecting edge endpoints.
206 blossom_children[blossom].clear();
207 blossom_endpoints[blossom].clear();
208 // Trace back from blossom_v to base.
209 while blossom_v != blossom_b {
210 // Add blossom_v to the new blossom
211 blossom_parents[blossom_v] = Some(blossom);
212 blossom_children[blossom].push(blossom_v);
213 let blossom_v_endpoint_label = label_ends[blossom_v].unwrap();
214 blossom_endpoints[blossom].push(blossom_v_endpoint_label);
215 assert!(
216 labels[blossom_v] == Some(2)
217 || (labels[blossom_v] == Some(1)
218 && label_ends[blossom_v]
219 == mate.get(&blossom_base[blossom_v].unwrap()).copied())
220 );
221 // Trace one step back.
222 assert!(label_ends[blossom_v].is_some());
223 v = endpoints[blossom_v_endpoint_label];
224 blossom_v = in_blossoms[v];
225 }
226 // Reverse lists, add endpoint that connects the pair of S vertices.
227 blossom_children[blossom].push(blossom_b);
228 blossom_children[blossom].reverse();
229 blossom_endpoints[blossom].reverse();
230 blossom_endpoints[blossom].push(2 * edge);
231 // Trace back from w to base.
232 while blossom_w != blossom_b {
233 // Add blossom_w to the new blossom
234 blossom_parents[blossom_w] = Some(blossom);
235 blossom_children[blossom].push(blossom_w);
236 let blossom_w_endpoint_label = label_ends[blossom_w].unwrap();
237 blossom_endpoints[blossom].push(blossom_w_endpoint_label ^ 1);
238 assert!(
239 labels[blossom_w] == Some(2)
240 || (labels[blossom_w] == Some(1)
241 && label_ends[blossom_w]
242 == mate.get(&blossom_base[blossom_w].unwrap()).copied())
243 );
244 // Trace one step back
245 assert!(label_ends[blossom_w].is_some());
246 w = endpoints[blossom_w_endpoint_label];
247 blossom_w = in_blossoms[w];
248 }
249 // Set label to S.
250 assert!(labels[blossom_b] == Some(1));
251 labels[blossom] = Some(1);
252 label_ends[blossom] = label_ends[blossom_b];
253 // Set dual variable to 0
254 dual_var[blossom] = 0;
255 // Relabel vertices
256 for node in blossom_leaves(blossom, num_nodes, blossom_children)? {
257 if labels[in_blossoms[node]] == Some(2) {
258 // This T-vertex now turns into an S-vertex because it becomes
259 // part of an S-blossom; add it to the queue
260 queue.push(node);
261 }
262 in_blossoms[node] = blossom;
263 }
264 // Compute blossom_best_edges[blossom]. Use IndexMap for deterministic iteration order.
265 let mut best_edge_to = DictMap::with_capacity(2 * num_nodes);
266 for bv_ref in &blossom_children[blossom] {
267 let bv = *bv_ref;
268 // This sub-blossom does not have a list of least-slack edges;
269 // get the information from the vertices.
270 let nblists: Vec<Vec<usize>> = if blossom_best_edges[bv].is_empty() {
271 let mut tmp: Vec<Vec<usize>> = Vec::new();
272 for node in blossom_leaves(bv, num_nodes, blossom_children)? {
273 tmp.push(neighbor_endpoints[node].iter().map(|p| p / 2).collect());
274 }
275 tmp
276 } else {
277 // Walk this sub-blossom's least-slack edges.
278 vec![blossom_best_edges[bv].clone()]
279 };
280 for nblist in nblists {
281 for edge_index in nblist {
282 let (mut i, mut j, _edge_weight) = edges[edge_index];
283 if in_blossoms[j] == blossom {
284 mem::swap(&mut i, &mut j);
285 }
286 let blossom_j = in_blossoms[j];
287 if blossom_j != blossom
288 && labels[blossom_j] == Some(1)
289 && (best_edge_to.get(&blossom_j).is_none()
290 || slack(edge_index, dual_var, edges)
291 < slack(best_edge_to[&blossom_j], dual_var, edges))
292 {
293 best_edge_to.insert(blossom_j, edge_index);
294 }
295 }
296 }
297 // Forget about least-slack edges of the sub-blossom.
298 blossom_best_edges[bv].clear();
299 best_edge[bv] = None;
300 }
301 blossom_best_edges[blossom] = best_edge_to.values().copied().collect();
302 //select best_edge[blossom]
303 best_edge[blossom] = None;
304 for edge_index in &blossom_best_edges[blossom] {
305 if best_edge[blossom].is_none()
306 || slack(*edge_index, dual_var, edges)
307 < slack(best_edge[blossom].unwrap(), dual_var, edges)
308 {
309 best_edge[blossom] = Some(*edge_index);
310 }
311 }
312 Ok(())
313}
314
315/// Expand the given top level blossom
316fn expand_blossom<E>(
317 blossom: usize,
318 end_stage: bool,
319 num_nodes: usize,
320 blossom_children: &mut Vec<Vec<usize>>,
321 blossom_parents: &mut Vec<Option<usize>>,
322 in_blossoms: &mut Vec<usize>,
323 dual_var: &[i128],
324 labels: &mut Vec<Option<usize>>,
325 label_ends: &mut Vec<Option<usize>>,
326 best_edge: &mut Vec<Option<usize>>,
327 queue: &mut Vec<usize>,
328 blossom_base: &mut Vec<Option<usize>>,
329 endpoints: &[usize],
330 mate: &HashMap<usize, usize>,
331 blossom_endpoints: &mut Vec<Vec<usize>>,
332 allowed_edge: &mut Vec<bool>,
333 unused_blossoms: &mut Vec<usize>,
334) -> Result<(), E> {
335 // Convert sub-blossoms into top-level blossoms.
336 for s in blossom_children[blossom].clone() {
337 blossom_parents[s] = None;
338 if s < num_nodes {
339 in_blossoms[s] = s
340 } else if end_stage && dual_var[s] == 0 {
341 // Recursively expand this sub-blossom
342 expand_blossom(
343 s,
344 end_stage,
345 num_nodes,
346 blossom_children,
347 blossom_parents,
348 in_blossoms,
349 dual_var,
350 labels,
351 label_ends,
352 best_edge,
353 queue,
354 blossom_base,
355 endpoints,
356 mate,
357 blossom_endpoints,
358 allowed_edge,
359 unused_blossoms,
360 )?;
361 } else {
362 for v in blossom_leaves(s, num_nodes, blossom_children)? {
363 in_blossoms[v] = s;
364 }
365 }
366 }
367 // if we expand a T-blossom during a stage, its a sub-blossoms must be
368 // relabeled
369 if !end_stage && labels[blossom] == Some(2) {
370 // start at the sub-blossom through which the expanding blossom
371 // obtained its label, and relabel sub-blossoms until we reach the
372 // base.
373 assert!(label_ends[blossom].is_some());
374 let entry_child = in_blossoms[endpoints[label_ends[blossom].unwrap() ^ 1]];
375 // Decide in which direction we will go around the blossom.
376 let i = blossom_children[blossom]
377 .iter()
378 .position(|x| *x == entry_child)
379 .unwrap();
380 let mut j = i as i128;
381 let j_step: i128;
382 let endpoint_trick: usize = if i & 1 != 0 {
383 // Start index is odd; go forward and wrap
384 j -= blossom_children[blossom].len() as i128;
385 j_step = 1;
386 0
387 } else {
388 // start index is even; go backward
389 j_step = -1;
390 1
391 };
392 // Move along the blossom until we get to the base
393 let mut p = label_ends[blossom].unwrap();
394 while j != 0 {
395 // Relabel the T-sub-blossom.
396 labels[endpoints[p ^ 1]] = Some(0);
397 if j < 0 {
398 let length = blossom_endpoints[blossom].len();
399 let index = length - j.unsigned_abs() as usize;
400 labels[endpoints
401 [blossom_endpoints[blossom][index - endpoint_trick] ^ endpoint_trick ^ 1]] =
402 Some(0);
403 } else {
404 labels[endpoints[blossom_endpoints[blossom][j as usize - endpoint_trick]
405 ^ endpoint_trick
406 ^ 1]] = Some(0);
407 }
408 assign_label(
409 endpoints[p ^ 1],
410 2,
411 Some(p),
412 num_nodes,
413 in_blossoms,
414 labels,
415 label_ends,
416 best_edge,
417 queue,
418 blossom_children,
419 blossom_base,
420 endpoints,
421 mate,
422 )?;
423 // Step to the next S-sub-blossom and note it's forward endpoint.
424 let endpoint_index = if j < 0 {
425 let tmp = j - endpoint_trick as i128;
426 let length = blossom_endpoints[blossom].len();
427 let index = length - tmp.unsigned_abs() as usize;
428 blossom_endpoints[blossom][index]
429 } else {
430 blossom_endpoints[blossom][j as usize - endpoint_trick]
431 };
432 allowed_edge[endpoint_index / 2] = true;
433 j += j_step;
434 p = if j < 0 {
435 let tmp = j - endpoint_trick as i128;
436 let length = blossom_endpoints[blossom].len();
437 let index = length - tmp.unsigned_abs() as usize;
438 blossom_endpoints[blossom][index] ^ endpoint_trick
439 } else {
440 blossom_endpoints[blossom][j as usize - endpoint_trick] ^ endpoint_trick
441 };
442 // Step to the next T-sub-blossom.
443 allowed_edge[p / 2] = true;
444 j += j_step;
445 }
446 // Relabel the base T-sub-blossom WITHOUT stepping through to
447 // its mate (so don't call assign_label())
448 let blossom_v = if j < 0 {
449 let length = blossom_children[blossom].len();
450 let index = length - j.unsigned_abs() as usize;
451 blossom_children[blossom][index]
452 } else {
453 blossom_children[blossom][j as usize]
454 };
455 labels[endpoints[p ^ 1]] = Some(2);
456 labels[blossom_v] = Some(2);
457 label_ends[endpoints[p ^ 1]] = Some(p);
458 label_ends[blossom_v] = Some(p);
459 best_edge[blossom_v] = None;
460 // Continue along the blossom until we get back to entry_child
461 j += j_step;
462 let mut j_index = if j < 0 {
463 let length = blossom_children[blossom].len();
464 length - j.unsigned_abs() as usize
465 } else {
466 j as usize
467 };
468 while blossom_children[blossom][j_index] != entry_child {
469 // Examine the vertices of the sub-blossom to see whether
470 // it is reachable from a neighboring S-vertex outside the
471 // expanding blososm.
472 let bv = blossom_children[blossom][j_index];
473 if labels[bv] == Some(1) {
474 // This sub-blossom just got label S through one of its
475 // neighbors; leave it
476 j += j_step;
477 if j < 0 {
478 let length = blossom_children[blossom].len();
479 j_index = length - j.unsigned_abs() as usize;
480 } else {
481 j_index = j as usize;
482 }
483 continue;
484 }
485 let mut v: usize = 0;
486 for tmp in blossom_leaves(bv, num_nodes, blossom_children)? {
487 v = tmp;
488 if labels[v].unwrap() != 0 {
489 break;
490 }
491 }
492 // If the sub-blossom contains a reachable vertex, assign label T
493 // to the sub-blossom
494 if labels[v] != Some(0) {
495 assert!(labels[v] == Some(2));
496 assert!(in_blossoms[v] == bv);
497 labels[v] = Some(0);
498 labels[endpoints[mate[&blossom_base[bv].unwrap()]]] = Some(0);
499 assign_label(
500 v,
501 2,
502 label_ends[v],
503 num_nodes,
504 in_blossoms,
505 labels,
506 label_ends,
507 best_edge,
508 queue,
509 blossom_children,
510 blossom_base,
511 endpoints,
512 mate,
513 )?;
514 }
515 j += j_step;
516 if j < 0 {
517 let length = blossom_children[blossom].len();
518 j_index = length - j.unsigned_abs() as usize;
519 } else {
520 j_index = j as usize;
521 }
522 }
523 }
524 // Recycle the blossom number.
525 labels[blossom] = None;
526 label_ends[blossom] = None;
527 blossom_children[blossom].clear();
528 blossom_endpoints[blossom].clear();
529 blossom_base[blossom] = None;
530 best_edge[blossom] = None;
531 unused_blossoms.push(blossom);
532 Ok(())
533}
534
535/// Swap matched/unmatched edges over an alternating path through blossom b
536/// between vertex v and the base vertex. Keep blossom bookkeeping consistent.
537fn augment_blossom(
538 blossom: usize,
539 node: usize,
540 num_nodes: usize,
541 blossom_parents: &[Option<usize>],
542 endpoints: &[usize],
543 blossom_children: &mut Vec<Vec<usize>>,
544 blossom_endpoints: &mut Vec<Vec<usize>>,
545 blossom_base: &mut Vec<Option<usize>>,
546 mate: &mut HashMap<usize, usize>,
547) {
548 // Bubble up through the blossom tree from vertex v to an immediate
549 // sub-blossom of b.
550 let mut tmp = node;
551 while blossom_parents[tmp] != Some(blossom) {
552 tmp = blossom_parents[tmp].unwrap();
553 }
554 // Recursively deal with the first sub-blossom.
555 if tmp >= num_nodes {
556 augment_blossom(
557 tmp,
558 node,
559 num_nodes,
560 blossom_parents,
561 endpoints,
562 blossom_children,
563 blossom_endpoints,
564 blossom_base,
565 mate,
566 );
567 }
568 // Decide in which direction we will go around the blossom.
569 let i = blossom_children[blossom]
570 .iter()
571 .position(|x| *x == tmp)
572 .unwrap();
573 let mut j: i128 = i as i128;
574 let j_step: i128;
575 let endpoint_trick: usize = if i & 1 != 0 {
576 // start index is odd; go forward and wrap.
577 j -= blossom_children[blossom].len() as i128;
578 j_step = 1;
579 0
580 } else {
581 // Start index is even; go backward.
582 j_step = -1;
583 1
584 };
585 // Move along the blossom until we get to the base.
586 while j != 0 {
587 // Step to the next sub-blossom and augment it recursively.
588 j += j_step;
589
590 tmp = if j < 0 {
591 let length = blossom_children[blossom].len();
592 let index = length - j.unsigned_abs() as usize;
593 blossom_children[blossom][index]
594 } else {
595 blossom_children[blossom][j as usize]
596 };
597 let p = if j < 0 {
598 let length = blossom_endpoints[blossom].len();
599 let index = length - j.unsigned_abs() as usize - endpoint_trick;
600 blossom_endpoints[blossom][index] ^ endpoint_trick
601 } else {
602 blossom_endpoints[blossom][j as usize - endpoint_trick] ^ endpoint_trick
603 };
604 if tmp > num_nodes {
605 augment_blossom(
606 tmp,
607 endpoints[p],
608 num_nodes,
609 blossom_parents,
610 endpoints,
611 blossom_children,
612 blossom_endpoints,
613 blossom_base,
614 mate,
615 );
616 }
617 j += j_step;
618 if j < 0 {
619 let length = blossom_children[blossom].len();
620 let index = length - j.unsigned_abs() as usize;
621 tmp = blossom_children[blossom][index];
622 } else {
623 tmp = blossom_children[blossom][j as usize];
624 }
625 if tmp > num_nodes {
626 augment_blossom(
627 tmp,
628 endpoints[p ^ 1],
629 num_nodes,
630 blossom_parents,
631 endpoints,
632 blossom_children,
633 blossom_endpoints,
634 blossom_base,
635 mate,
636 );
637 }
638 // Match the edge connecting those sub-blossoms.
639 mate.insert(endpoints[p], p ^ 1);
640 mate.insert(endpoints[p ^ 1], p);
641 }
642 // Rotate the list of sub-blossoms to put the new base at the front.
643 let mut children: Vec<usize> = blossom_children[blossom][i..].to_vec();
644 children.extend_from_slice(&blossom_children[blossom][..i]);
645 blossom_children[blossom] = children;
646 let mut endpoint: Vec<usize> = blossom_endpoints[blossom][i..].to_vec();
647 endpoint.extend_from_slice(&blossom_endpoints[blossom][..i]);
648 blossom_endpoints[blossom] = endpoint;
649 blossom_base[blossom] = blossom_base[blossom_children[blossom][0]];
650 assert!(blossom_base[blossom] == Some(node));
651}
652
653/// Swap matched/unmatched edges over an alternating path between two
654/// single vertices. The augmenting path runs through edge k, which
655/// connects a pair of S vertices.
656fn augment_matching(
657 edge: usize,
658 num_nodes: usize,
659 edges: &[(usize, usize, i128)],
660 in_blossoms: &[usize],
661 labels: &[Option<usize>],
662 label_ends: &[Option<usize>],
663 blossom_parents: &[Option<usize>],
664 endpoints: &[usize],
665 blossom_children: &mut Vec<Vec<usize>>,
666 blossom_endpoints: &mut Vec<Vec<usize>>,
667 blossom_base: &mut Vec<Option<usize>>,
668 mate: &mut HashMap<usize, usize>,
669) {
670 let (v, w, _weight) = edges[edge];
671 for (s_ref, p_ref) in [(v, 2 * edge + 1), (w, 2 * edge)].iter() {
672 // Match vertex s to remote endpoint p. Then trace back from s
673 // until we find a single vertex, swapping matched and unmatched
674 // edges as we go.
675 let mut s: usize = *s_ref;
676 let mut p: usize = *p_ref;
677 loop {
678 let blossom_s = in_blossoms[s];
679 assert!(labels[blossom_s] == Some(1));
680 assert!(label_ends[blossom_s] == mate.get(&blossom_base[blossom_s].unwrap()).copied());
681 // Augment through the S-blossom from s to base.
682 if blossom_s >= num_nodes {
683 augment_blossom(
684 blossom_s,
685 s,
686 num_nodes,
687 blossom_parents,
688 endpoints,
689 blossom_children,
690 blossom_endpoints,
691 blossom_base,
692 mate,
693 );
694 }
695 // Update mate[s]
696 mate.insert(s, p);
697 // Trace one step back.
698 if label_ends[blossom_s].is_none() {
699 // Reached a single vertex; stop
700 break;
701 }
702 let t = endpoints[label_ends[blossom_s].unwrap()];
703 let blossom_t = in_blossoms[t];
704 assert!(labels[blossom_t] == Some(2));
705 // Trace one step back
706 assert!(label_ends[blossom_t].is_some());
707 s = endpoints[label_ends[blossom_t].unwrap()];
708 let j = endpoints[label_ends[blossom_t].unwrap() ^ 1];
709 // Augment through the T-blossom from j to base.
710 assert!(blossom_base[blossom_t] == Some(t));
711 if blossom_t >= num_nodes {
712 augment_blossom(
713 blossom_t,
714 j,
715 num_nodes,
716 blossom_parents,
717 endpoints,
718 blossom_children,
719 blossom_endpoints,
720 blossom_base,
721 mate,
722 );
723 }
724 // Update mate[j]
725 mate.insert(j, label_ends[blossom_t].unwrap());
726 // Keep the opposite endpoint;
727 // it will be assigned to mate[s] in the next step.
728 p = label_ends[blossom_t].unwrap() ^ 1;
729 }
730 }
731}
732
733/// Swap matched/unmatched edges over an alternating path between two
734/// single vertices. The augmenting path runs through the edge, which
735/// connects a pair of S vertices.
736fn verify_optimum(
737 max_cardinality: bool,
738 num_nodes: usize,
739 num_edges: usize,
740 edges: &[(usize, usize, i128)],
741 endpoints: &[usize],
742 dual_var: &[i128],
743 blossom_parents: &[Option<usize>],
744 blossom_endpoints: &[Vec<usize>],
745 blossom_base: &[Option<usize>],
746 mate: &HashMap<usize, usize>,
747) {
748 let dual_var_node_min: i128 = *dual_var[..num_nodes].iter().min().unwrap();
749 let node_dual_offset: i128 = if max_cardinality {
750 // Vertices may have negative dual;
751 // find a constant non-negative number to add to all vertex duals.
752 max(0, -dual_var_node_min)
753 } else {
754 0
755 };
756 assert!(dual_var_node_min + node_dual_offset >= 0);
757 assert!(*dual_var[num_nodes..].iter().min().unwrap() >= 0);
758 // 0. all edges have non-negative slack and
759 // 1. all matched edges have zero slack;
760 for (edge, (i, j, weight)) in edges.iter().enumerate().take(num_edges) {
761 let mut s = dual_var[*i] + dual_var[*j] - 2 * weight;
762 let mut i_blossoms: Vec<usize> = vec![*i];
763 let mut j_blossoms: Vec<usize> = vec![*j];
764 while blossom_parents[*i_blossoms.last().unwrap()].is_some() {
765 i_blossoms.push(blossom_parents[*i_blossoms.last().unwrap()].unwrap());
766 }
767 while blossom_parents[*j_blossoms.last().unwrap()].is_some() {
768 j_blossoms.push(blossom_parents[*j_blossoms.last().unwrap()].unwrap());
769 }
770 i_blossoms.reverse();
771 j_blossoms.reverse();
772 for (blossom_i, blossom_j) in i_blossoms.iter().zip(j_blossoms.iter()) {
773 if blossom_i != blossom_j {
774 break;
775 }
776 s += 2 * dual_var[*blossom_i];
777 }
778 assert!(s >= 0);
779
780 if (mate.get(i).is_some() && mate.get(i).unwrap() / 2 == edge)
781 || (mate.get(j).is_some() && mate.get(j).unwrap() / 2 == edge)
782 {
783 assert!(mate[i] / 2 == edge && mate[j] / 2 == edge);
784 assert!(s == 0);
785 }
786 }
787 // 2. all single vertices have zero dual value;
788 for (node, dual_var_node) in dual_var.iter().enumerate().take(num_nodes) {
789 assert!(mate.get(&node).is_some() || dual_var_node + node_dual_offset == 0);
790 }
791 // 3. all blossoms with positive dual value are full.
792 for blossom in num_nodes..2 * num_nodes {
793 if blossom_base[blossom].is_some() && dual_var[blossom] > 0 {
794 assert!(blossom_endpoints[blossom].len() % 2 == 1);
795 for p in blossom_endpoints[blossom].iter().skip(1).step_by(2) {
796 assert!(mate.get(&endpoints[*p]).copied() == Some(p ^ 1));
797 assert!(mate.get(&endpoints[*p ^ 1]) == Some(p));
798 }
799 }
800 }
801}
802
803/// Compute a maximum-weighted matching in the general undirected weighted
804/// graph given by "edges". If `max_cardinality` is true, only
805/// maximum-cardinality matchings are considered as solutions.
806///
807/// The algorithm is based on "Efficient Algorithms for Finding Maximum
808/// Matching in Graphs" by Zvi Galil, ACM Computing Surveys, 1986.
809///
810/// Based on networkx implementation
811/// <https://github.com/networkx/networkx/blob/3351206a3ce5b3a39bb2fc451e93ef545b96c95b/networkx/algorithms/matching.py>
812///
813/// With reference to the standalone prototype implementation from:
814/// <http://jorisvr.nl/article/maximum-matching>
815///
816/// <http://jorisvr.nl/files/graphmatching/20130407/mwmatching.py>
817///
818/// The function takes time O(n**3)
819///
820/// Arguments:
821///
822/// * `graph` - The undirected graph to compute the maximum weight matching for
823/// * `max_cardinality` - If set to true compute the maximum-cardinality matching
824/// with maximum weight among all maximum-cardinality matchings
825/// * `weight_fn` - A callback function that will be give a edge reference and
826/// expected to return a `i128` representing the weight of the edge
827/// * `verify_optimum_flag`: If true an prior to returning an additional routine
828/// to verify the optimal solution was found will be run after computing
829/// the maximum weight matching. If it's true and the found matching is not
830/// an optimal solution this function will panic. This option should
831/// normally be only set true during testing.
832///
833/// # Example
834/// ```rust
835/// use rustworkx_core::petgraph;
836/// use rustworkx_core::max_weight_matching::max_weight_matching;
837/// use rustworkx_core::Result;
838///
839/// use hashbrown::HashSet;
840///
841/// // Create a path graph
842/// let g = petgraph::graph::UnGraph::<i32, i128>::from_edges(&[
843/// (1, 2, 5), (2, 3, 11), (3, 4, 5)
844/// ]);
845///
846/// // Run max weight matching with max cardinality set to false
847/// let res: Result<HashSet<(usize, usize)>> = max_weight_matching(
848/// &g, false, |e| Ok(*e.weight()), true
849/// );
850/// // Run max weight matching with max cardinality set to true
851/// let maxc_res: Result<HashSet<(usize, usize)>> = max_weight_matching(
852/// &g, true, |e| Ok(*e.weight()), true
853/// );
854///
855/// let matching = res.unwrap();
856/// let maxc_matching = maxc_res.unwrap();
857/// // Check output
858/// assert_eq!(matching.len(), 1);
859/// assert!(matching.contains(&(2, 3)) || matching.contains(&(3, 2)));
860/// assert_eq!(maxc_matching.len(), 2);
861/// assert!(maxc_matching.contains(&(1, 2)) || maxc_matching.contains(&(2, 1)));
862/// assert!(maxc_matching.contains(&(3, 4)) || maxc_matching.contains(&(4, 3)));
863/// ```
864pub fn max_weight_matching<G, F, E>(
865 graph: G,
866 max_cardinality: bool,
867 mut weight_fn: F,
868 verify_optimum_flag: bool,
869) -> Result<HashSet<(usize, usize)>, E>
870where
871 G: EdgeCount
872 + NodeCount
873 + IntoNodeIdentifiers
874 + GraphProp<EdgeType = Undirected>
875 + GraphBase
876 + NodeIndexable
877 + IntoEdges,
878 F: FnMut(G::EdgeRef) -> Result<i128, E>,
879 <G as GraphBase>::NodeId: std::cmp::Eq + Hash,
880{
881 let num_edges = graph.edge_count();
882 let num_nodes = graph.node_count();
883 let mut out_set: HashSet<(usize, usize)> = HashSet::with_capacity(num_nodes);
884 // Exit fast for graph without edges
885 if num_edges == 0 {
886 return Ok(HashSet::new());
887 }
888 // Node indices in the PyGraph may not be contiguous however the
889 // algorithm operates on contiguous indices 0..num_nodes. node_map maps
890 // the PyGraph's NodeIndex to the contiguous usize used inside the
891 // algorithm
892 let node_map: HashMap<G::NodeId, usize> = graph
893 .node_identifiers()
894 .enumerate()
895 .map(|(index, node_index)| (node_index, index))
896 .collect();
897 let mut edges: Vec<(usize, usize, i128)> = Vec::with_capacity(num_edges);
898 let mut max_weight: i128 = 0;
899 for edge in graph.edge_references() {
900 let edge_weight: i128 = weight_fn(edge)?;
901 if edge_weight > max_weight {
902 max_weight = edge_weight;
903 };
904 edges.push((
905 node_map[&edge.source()],
906 node_map[&edge.target()],
907 edge_weight,
908 ));
909 }
910 // If p is an edge endpoint
911 // endpoints[p] is the node index to which endpoint p is attached
912 let endpoints: Vec<usize> = (0..2 * num_edges)
913 .map(|endpoint| {
914 let edge_tuple = edges[endpoint / 2];
915 let out_value: usize = if endpoint % 2 == 0 {
916 edge_tuple.0
917 } else {
918 edge_tuple.1
919 };
920 out_value
921 })
922 .collect();
923 // If v is a node/vertex
924 // neighbor_endpoints[v] is the list of remote endpoints of the edges
925 // attached to v.
926 // Not modified by algorithm (only mut to initially construct contents).
927 let mut neighbor_endpoints: Vec<Vec<usize>> = (0..num_nodes).map(|_| Vec::new()).collect();
928 for edge in 0..num_edges {
929 neighbor_endpoints[edges[edge].0].push(2 * edge + 1);
930 neighbor_endpoints[edges[edge].1].push(2 * edge);
931 }
932 // If v is a vertex,
933 // mate[v] is the remote endpoint of its matched edge, or None if it is
934 // single (i.e. endpoint[mate[v]] is v's partner vertex).
935 // Initially all vertices are single; updated during augmentation.
936 let mut mate: HashMap<usize, usize> = HashMap::with_capacity(num_nodes);
937 // If b is a top-level blossom
938 // label[b] is 0 if b is unlabeled (free);
939 // 1 if b is a S-vertex/blossom;
940 // 2 if b is a T-vertex/blossom;
941 // The label of a vertex/node is found by looking at the label of its
942 // top-level containing blossom
943 // If v is a node/vertex inside a T-Blossom,
944 // label[v] is 2 if and only if v is reachable from an S_Vertex outside
945 // the blossom.
946 let mut labels: Vec<Option<usize>>;
947 // If b is a labeled top-level blossom,
948 // label_ends[b] is the remote endpoint of the edge through which b
949 // obtained its label, or None if b's base vertex is single.
950 // If v is a vertex inside a T-blossom and label[v] == 2,
951 // label_ends[v] is the remote endpoint of the edge through which v is
952 // reachable from outside the blossom
953 let mut label_ends: Vec<Option<usize>>;
954 // If v is a vertex/node
955 // in_blossoms[v] is the top-level blossom to which v belongs.
956 // If v is a top-level vertex, v is itself a blossom (a trivial blossom)
957 // and in_blossoms[v] == v.
958 // Initially all nodes are top-level trivial blossoms.
959 let mut in_blossoms: Vec<usize> = (0..num_nodes).collect();
960 // if b is a sub-blossom
961 // blossom_parents[b] is its immediate parent
962 // If b is a top-level blossom, blossom_parents[b] is None
963 let mut blossom_parents: Vec<Option<usize>> = vec![None; 2 * num_nodes];
964 // If b is a non-trivial (sub-)blossom,
965 // blossom_children[b] is an ordered list of its sub-blossoms, starting with
966 // the base and going round the blossom.
967 let mut blossom_children: Vec<Vec<usize>> = (0..2 * num_nodes).map(|_| Vec::new()).collect();
968 // If b is a (sub-)blossom,
969 // blossombase[b] is its base VERTEX (i.e. recursive sub-blossom).
970 let mut blossom_base: Vec<Option<usize>> = (0..num_nodes).map(Some).collect();
971 blossom_base.append(&mut vec![None; num_nodes]);
972 // If b is a non-trivial (sub-)blossom,
973 // blossom_endpoints[b] is a list of endpoints on its connecting edges,
974 // such that blossom_endpoints[b][i] is the local endpoint of
975 // blossom_children[b][i] on the edge that connects it to
976 // blossom_children[b][wrap(i+1)].
977 let mut blossom_endpoints: Vec<Vec<usize>> = (0..2 * num_nodes).map(|_| Vec::new()).collect();
978 // If v is a free vertex (or an unreached vertex inside a T-blossom),
979 // best_edge[v] is the edge to an S-vertex with least slack,
980 // or None if there is no such edge. If b is a (possibly trivial)
981 // top-level S-blossom,
982 // best_edge[b] is the least-slack edge to a different S-blossom,
983 // or None if there is no such edge.
984 let mut best_edge: Vec<Option<usize>>;
985 // If b is a non-trivial top-level S-blossom,
986 // blossom_best_edges[b] is a list of least-slack edges to neighboring
987 // S-blossoms, or None if no such list has been computed yet.
988 // This is used for efficient computation of delta3.
989 let mut blossom_best_edges: Vec<Vec<usize>> = (0..2 * num_nodes).map(|_| Vec::new()).collect();
990 let mut unused_blossoms: Vec<usize> = (num_nodes..2 * num_nodes).collect();
991 // If v is a vertex,
992 // dual_var[v] = 2 * u(v) where u(v) is the v's variable in the dual
993 // optimization problem (multiplication by two ensures integer values
994 // throughout the algorithm if all edge weights are integers).
995 // If b is a non-trivial blossom,
996 // dual_var[b] = z(b) where z(b) is b's variable in the dual optimization
997 // problem.
998 // dual_var is for vertex v in 0..num_nodes and blossom b in
999 // num_nodes..2* num nodes
1000 let mut dual_var: Vec<i128> = vec![max_weight; num_nodes];
1001 dual_var.append(&mut vec![0; num_nodes]);
1002 // If allowed_edge[k] is true, edge k has zero slack in the optimization
1003 // problem; if allowed_edge[k] is false, the edge's slack may or may not
1004 // be zero.
1005 let mut allowed_edge: Vec<bool>;
1006 // Queue of newly discovered S-vertices
1007 let mut queue: Vec<usize> = Vec::with_capacity(num_nodes);
1008
1009 // Main loop: continue until no further improvement is possible
1010 for _ in 0..num_nodes {
1011 // Each iteration of this loop is a "stage".
1012 // A stage finds an augmenting path and uses that to improve
1013 // the matching.
1014
1015 // Removal labels from top-level blossoms/vertices
1016 labels = vec![Some(0); 2 * num_nodes];
1017 label_ends = vec![None; 2 * num_nodes];
1018 // Forget all about least-slack edges.
1019 best_edge = vec![None; 2 * num_nodes];
1020 blossom_best_edges.splice(num_nodes.., (0..num_nodes).map(|_| Vec::new()));
1021 // Loss of labeling means that we can not be sure that currently
1022 // allowable edges remain allowable throughout this stage.
1023 allowed_edge = vec![false; num_edges];
1024 // Make queue empty
1025 queue.clear();
1026 // Label single blossom/vertices with S and put them in queue.
1027 for v in 0..num_nodes {
1028 if mate.get(&v).is_none() && labels[in_blossoms[v]] == Some(0) {
1029 assign_label(
1030 v,
1031 1,
1032 None,
1033 num_nodes,
1034 &in_blossoms,
1035 &mut labels,
1036 &mut label_ends,
1037 &mut best_edge,
1038 &mut queue,
1039 &blossom_children,
1040 &blossom_base,
1041 &endpoints,
1042 &mate,
1043 )?;
1044 }
1045 }
1046 // Loop until we succeed in augmenting the matching.
1047 let mut augmented = false;
1048 loop {
1049 // Each iteration of this loop is a "substage".
1050 // A substage tries to find an augmenting path;
1051 // if found, the path is used to improve the matching and
1052 // the stage ends. If there is no augmenting path, the
1053 // primal-dual method is used to find some slack out of
1054 // the dual variables.
1055
1056 // Continue labeling until all vertices which are reachable
1057 // through an alternating path have got a label.
1058 while !queue.is_empty() && !augmented {
1059 // Take an S vertex from the queue
1060 let v = queue.pop().unwrap();
1061 assert!(labels[in_blossoms[v]] == Some(1));
1062
1063 // Scan its neighbors
1064 for p in &neighbor_endpoints[v] {
1065 let k = *p / 2;
1066 let mut kslack = 0;
1067 let w = endpoints[*p];
1068 // w is a neighbor of v
1069 if in_blossoms[v] == in_blossoms[w] {
1070 // this edge is internal to a blossom; ignore it
1071 continue;
1072 }
1073 if !allowed_edge[k] {
1074 kslack = slack(k, &dual_var, &edges);
1075 if kslack <= 0 {
1076 // edge k has zero slack -> it is allowable
1077 allowed_edge[k] = true;
1078 }
1079 }
1080 if allowed_edge[k] {
1081 if labels[in_blossoms[w]] == Some(0) {
1082 // (C1) w is a free vertex;
1083 // label w with T and label its mate with S (R12).
1084 assign_label(
1085 w,
1086 2,
1087 Some(*p ^ 1),
1088 num_nodes,
1089 &in_blossoms,
1090 &mut labels,
1091 &mut label_ends,
1092 &mut best_edge,
1093 &mut queue,
1094 &blossom_children,
1095 &blossom_base,
1096 &endpoints,
1097 &mate,
1098 )?;
1099 } else if labels[in_blossoms[w]] == Some(1) {
1100 // (C2) w is an S-vertex (not in the same blossom);
1101 // follow back-links to discover either an
1102 // augmenting path or a new blossom.
1103 let base = scan_blossom(
1104 v,
1105 w,
1106 &in_blossoms,
1107 &blossom_base,
1108 &endpoints,
1109 &mut labels,
1110 &label_ends,
1111 &mate,
1112 );
1113 match base {
1114 // Found a new blossom; add it to the blossom
1115 // bookkeeping and turn it into an S-blossom.
1116 Some(base) => add_blossom(
1117 base,
1118 k,
1119 &mut blossom_children,
1120 num_nodes,
1121 &edges,
1122 &mut in_blossoms,
1123 &mut dual_var,
1124 &mut labels,
1125 &mut label_ends,
1126 &mut best_edge,
1127 &mut queue,
1128 &mut blossom_base,
1129 &endpoints,
1130 &mut blossom_endpoints,
1131 &mut unused_blossoms,
1132 &mut blossom_best_edges,
1133 &mut blossom_parents,
1134 &neighbor_endpoints,
1135 &mate,
1136 )?,
1137 // Found an augmenting path; augment the
1138 // matching and end this stage.
1139 None => {
1140 augment_matching(
1141 k,
1142 num_nodes,
1143 &edges,
1144 &in_blossoms,
1145 &labels,
1146 &label_ends,
1147 &blossom_parents,
1148 &endpoints,
1149 &mut blossom_children,
1150 &mut blossom_endpoints,
1151 &mut blossom_base,
1152 &mut mate,
1153 );
1154 augmented = true;
1155 break;
1156 }
1157 };
1158 } else if labels[w] == Some(0) {
1159 // w is inside a T-blossom, but w itself has not
1160 // yet been reached from outside the blossom;
1161 // mark it as reached (we need this to relabel
1162 // during T-blossom expansion).
1163 assert!(labels[in_blossoms[w]] == Some(2));
1164 labels[w] = Some(2);
1165 label_ends[w] = Some(*p ^ 1);
1166 }
1167 } else if labels[in_blossoms[w]] == Some(1) {
1168 // Keep track of the least-slack non-allowable edge to
1169 // a different S-blossom
1170 let blossom = in_blossoms[v];
1171 if best_edge[blossom].is_none()
1172 || kslack < slack(best_edge[blossom].unwrap(), &dual_var, &edges)
1173 {
1174 best_edge[blossom] = Some(k);
1175 }
1176 } else if labels[w] == Some(0) {
1177 // w is a free vertex (or an unreached vertex inside
1178 // a T-blossom) but we can not reach it yet;
1179 // keep track of the least-slack edge that reaches w.
1180 if best_edge[w].is_none()
1181 || kslack < slack(best_edge[w].unwrap(), &dual_var, &edges)
1182 {
1183 best_edge[w] = Some(k)
1184 }
1185 }
1186 }
1187 }
1188 if augmented {
1189 break;
1190 }
1191 // There is no augmenting path under these constraints;
1192 // compute delta and reduce slack in the optimization problem.
1193 // (Note that our vertex dual variables, edge slacks and delta's
1194 // are pre-multiplied by two.)
1195 let mut delta_type = -1;
1196 let mut delta: Option<i128> = None;
1197 let mut delta_edge: Option<usize> = None;
1198 let mut delta_blossom: Option<usize> = None;
1199
1200 // Compute delta1: the minimum value of any vertex dual.
1201 if !max_cardinality {
1202 delta_type = 1;
1203 delta = Some(*dual_var[..num_nodes].iter().min().unwrap());
1204 }
1205
1206 // Compute delta2: the minimum slack on any edge between
1207 // an S-vertex and a free vertex.
1208 for v in 0..num_nodes {
1209 if labels[in_blossoms[v]] == Some(0) && best_edge[v].is_some() {
1210 let d = slack(best_edge[v].unwrap(), &dual_var, &edges);
1211 if delta_type == -1 || Some(d) < delta {
1212 delta = Some(d);
1213 delta_type = 2;
1214 delta_edge = best_edge[v];
1215 }
1216 }
1217 }
1218
1219 // Compute delta3: half the minimum slack on any edge between a
1220 // pair of S-blossoms
1221 for blossom in 0..2 * num_nodes {
1222 if blossom_parents[blossom].is_none()
1223 && labels[blossom] == Some(1)
1224 && best_edge[blossom].is_some()
1225 {
1226 let kslack = slack(best_edge[blossom].unwrap(), &dual_var, &edges);
1227 assert!(kslack % 2 == 0);
1228 let d = Some(kslack / 2);
1229 if delta_type == -1 || d < delta {
1230 delta = d;
1231 delta_type = 3;
1232 delta_edge = best_edge[blossom];
1233 }
1234 }
1235 }
1236
1237 // Compute delta4: minimum z variable of any T-blossom
1238 for blossom in num_nodes..2 * num_nodes {
1239 if blossom_base[blossom].is_some()
1240 && blossom_parents[blossom].is_none()
1241 && labels[blossom] == Some(2)
1242 && (delta_type == -1 || dual_var[blossom] < delta.unwrap())
1243 {
1244 delta = Some(dual_var[blossom]);
1245 delta_type = 4;
1246 delta_blossom = Some(blossom);
1247 }
1248 }
1249 if delta_type == -1 {
1250 // No further improvement possible; max-cardinality optimum
1251 // reached. Do a final delta update to make the optimum
1252 // verifiable
1253 assert!(max_cardinality);
1254 delta_type = 1;
1255 delta = Some(max(0, *dual_var[..num_nodes].iter().min().unwrap()));
1256 }
1257
1258 // Update dual variables according to delta.
1259 for v in 0..num_nodes {
1260 if labels[in_blossoms[v]] == Some(1) {
1261 // S-vertex: 2*u = 2*u - 2*delta
1262 dual_var[v] -= delta.unwrap();
1263 } else if labels[in_blossoms[v]] == Some(2) {
1264 // T-vertex: 2*u = 2*u + 2*delta
1265 dual_var[v] += delta.unwrap();
1266 }
1267 }
1268 for b in num_nodes..2 * num_nodes {
1269 if blossom_base[b].is_some() && blossom_parents[b].is_none() {
1270 if labels[b] == Some(1) {
1271 // top-level S-blossom: z = z + 2*delta
1272 dual_var[b] += delta.unwrap();
1273 } else if labels[b] == Some(2) {
1274 // top-level T-blossom: z = z - 2*delta
1275 dual_var[b] -= delta.unwrap();
1276 }
1277 }
1278 }
1279 // Take action at the point where minimum delta occurred.
1280 if delta_type == 1 {
1281 // No further improvement possible; optimum reached
1282 break;
1283 } else if delta_type == 2 {
1284 // Use the least-slack edge to continue the search.
1285 allowed_edge[delta_edge.unwrap()] = true;
1286 let (mut i, mut j, _weight) = edges[delta_edge.unwrap()];
1287 if labels[in_blossoms[i]] == Some(0) {
1288 mem::swap(&mut i, &mut j);
1289 }
1290 assert!(labels[in_blossoms[i]] == Some(1));
1291 queue.push(i);
1292 } else if delta_type == 3 {
1293 // Use the least-slack edge to continue the search.
1294 allowed_edge[delta_edge.unwrap()] = true;
1295 let (i, _j, _weight) = edges[delta_edge.unwrap()];
1296 assert!(labels[in_blossoms[i]] == Some(1));
1297 queue.push(i);
1298 } else if delta_type == 4 {
1299 // Expand the least-z blossom
1300 expand_blossom(
1301 delta_blossom.unwrap(),
1302 false,
1303 num_nodes,
1304 &mut blossom_children,
1305 &mut blossom_parents,
1306 &mut in_blossoms,
1307 &dual_var,
1308 &mut labels,
1309 &mut label_ends,
1310 &mut best_edge,
1311 &mut queue,
1312 &mut blossom_base,
1313 &endpoints,
1314 &mate,
1315 &mut blossom_endpoints,
1316 &mut allowed_edge,
1317 &mut unused_blossoms,
1318 )?;
1319 }
1320 // end of this substage
1321 }
1322 // Stop when no more augment paths can be found
1323 if !augmented {
1324 break;
1325 }
1326
1327 // End of a stage; expand all S-blossoms which have a dual_var == 0
1328 for blossom in num_nodes..2 * num_nodes {
1329 if blossom_parents[blossom].is_none()
1330 && blossom_base[blossom].is_some()
1331 && labels[blossom] == Some(1)
1332 && dual_var[blossom] == 0
1333 {
1334 expand_blossom(
1335 blossom,
1336 true,
1337 num_nodes,
1338 &mut blossom_children,
1339 &mut blossom_parents,
1340 &mut in_blossoms,
1341 &dual_var,
1342 &mut labels,
1343 &mut label_ends,
1344 &mut best_edge,
1345 &mut queue,
1346 &mut blossom_base,
1347 &endpoints,
1348 &mate,
1349 &mut blossom_endpoints,
1350 &mut allowed_edge,
1351 &mut unused_blossoms,
1352 )?;
1353 }
1354 }
1355 }
1356 if verify_optimum_flag {
1357 verify_optimum(
1358 max_cardinality,
1359 num_nodes,
1360 num_edges,
1361 &edges,
1362 &endpoints,
1363 &dual_var,
1364 &blossom_parents,
1365 &blossom_endpoints,
1366 &blossom_base,
1367 &mate,
1368 );
1369 }
1370
1371 // Transform mate[] such that mate[v] is the vertex to which v is paired
1372 // Also handle holes in node indices from PyGraph node removals by mapping
1373 // linear index to node index.
1374 let mut seen: HashSet<(usize, usize)> = HashSet::with_capacity(2 * num_nodes);
1375 let node_list: Vec<G::NodeId> = graph.node_identifiers().collect();
1376 for (index, node) in mate.iter() {
1377 let tmp = (
1378 graph.to_index(node_list[*index]),
1379 graph.to_index(node_list[endpoints[*node]]),
1380 );
1381 let rev_tmp = (
1382 graph.to_index(node_list[endpoints[*node]]),
1383 graph.to_index(node_list[*index]),
1384 );
1385 if !seen.contains(&tmp) && !seen.contains(&rev_tmp) {
1386 out_set.insert(tmp);
1387 seen.insert(tmp);
1388 seen.insert(rev_tmp);
1389 }
1390 }
1391 Ok(out_set)
1392}