From 98fa26eb2fbbcd416f1c8ca81f2010ca65629e61 Mon Sep 17 00:00:00 2001 From: Matthew Treinish Date: Fri, 15 Jan 2021 17:16:58 -0500 Subject: [PATCH] Add max weight matching function MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit This commit adds a new python function max_weight_function() for computing the maximum-weighted matching of a PyGraph object. The implementation of this function is based on “Efficient Algorithms for Finding Maximum Matching in Graphs”, Zvi Galil, ACM Computing Surveys, 1986. [1] It is basically a porting of the networkx implementation of the algorithm [2] with some additional inspiration from the prototype for the networkx version (it's basically the same code but some aspects were a bit clearer to figure out from the prototype rather than networkx's copy of it). [3][4] Fixes #216 [1] https://dl.acm.org/doi/10.1145/6462.6502 [2] https://github.com/networkx/networkx/blob/3351206a3ce5b3a39bb2fc451e93ef545b96c95b/networkx/algorithms/matching.py [3] http://jorisvr.nl/article/maximum-matching [4] http://jorisvr.nl/files/graphmatching/20130407/mwmatching.py --- Cargo.lock | 8 +- docs/source/api.rst | 1 + ...-max-weight-matching-089ba67668b3b071.yaml | 5 + src/lib.rs | 57 + src/max_weight_matching.rs | 1410 +++++++++++++++++ tests/test_max_weight_matching.py | 303 ++++ 6 files changed, 1780 insertions(+), 4 deletions(-) create mode 100644 releasenotes/notes/add-max-weight-matching-089ba67668b3b071.yaml create mode 100644 src/max_weight_matching.rs create mode 100644 tests/test_max_weight_matching.py diff --git a/Cargo.lock b/Cargo.lock index 91afda89c..a4c3302c1 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -78,9 +78,9 @@ dependencies = [ [[package]] name = "ctor" -version = "0.1.17" +version = "0.1.18" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "373c88d9506e2e9230f6107701b7d8425f4cb3f6df108ec3042a26e936666da5" +checksum = "10bcb9d7dcbf7002aaffbb53eac22906b64cdcc127971dcc387d8eb7c95d5560" dependencies = [ "quote", "syn", @@ -132,9 +132,9 @@ dependencies = [ [[package]] name = "hermit-abi" -version = "0.1.17" +version = "0.1.18" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "5aca5565f760fb5b220e499d72710ed156fdb74e631659e99377d9ebfbd13ae8" +checksum = "322f4de77956e22ed0e5032c359a0f1273f1f7f0d79bfa3b8ffbc730d7fbcc5c" dependencies = [ "libc", ] diff --git a/docs/source/api.rst b/docs/source/api.rst index 0ac5fcdb7..79694d06c 100644 --- a/docs/source/api.rst +++ b/docs/source/api.rst @@ -87,6 +87,7 @@ Algorithm Functions retworkx.digraph_dfs_edges retworkx.digraph_find_cycle retworkx.digraph_union + retworkx.max_weight_matching Exceptions ---------- diff --git a/releasenotes/notes/add-max-weight-matching-089ba67668b3b071.yaml b/releasenotes/notes/add-max-weight-matching-089ba67668b3b071.yaml new file mode 100644 index 000000000..0535ad677 --- /dev/null +++ b/releasenotes/notes/add-max-weight-matching-089ba67668b3b071.yaml @@ -0,0 +1,5 @@ +--- +features: + - | + Add a new function, :func:`~retworkx.max_weight_matching` for computing the + maximum-weighted matching for a :class:`~retworkx.PyGraph` object. diff --git a/src/lib.rs b/src/lib.rs index 9fb9b7fa5..28276ab4f 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -19,6 +19,7 @@ mod generators; mod graph; mod iterators; mod k_shortest_path; +mod max_weight_matching; mod union; use std::cmp::{Ordering, Reverse}; @@ -2475,6 +2476,61 @@ pub fn cycle_basis( cycles } +/// Compute a maximum-weighted matching for a :class:`~retworkx.PyGraph` +/// +/// A matching is a subset of edges in which no node occurs more than once. +/// The weight of a matching is the sum of the weights of its edges. +/// A maximal matching cannot add more edges and still be a matching. +/// The cardinality of a matching is the number of matched edges. +/// +/// This function takes time :math:`O(n^3)` where ``n` is the number of nodes +/// in the graph. +/// +/// This method is based on the "blossom" method for finding augmenting +/// paths and the "primal-dual" method for finding a matching of maximum +/// weight, both methods invented by Jack Edmonds [1]_. +/// +/// :param PyGraph graph: The undirected graph to compute the max weight +/// matching for. Expects to have no parallel edges (multigraphs are +/// untested currently). +/// :param max_cardinality bool: If True, compute the maximum-cardinality +/// matching with maximum weight among all maximum-cardinality matchings. +/// Defaults False. +/// :param callable weight_fn: An optional callable that will be passed a +/// single argument the edge object for each edge in the graph. It is +/// expected to return an ``int`` weight for that edge. For example, +/// if the weights are all integers you can use: ``lambda x: x``. If not +/// specified the value for ``default_weight`` will be used for all +/// edge weights. +/// :param int default_weight: The ``int`` value to use for all edge weights +/// in the graph if ``weight_fn`` is not specified. Defaults to ``1``. +/// +/// :returns: A dictionary of matching, the key is the node index and the +/// value is it's matched node index. Note that both directions will be +/// listed in the output, for example: ``{0: 1, 1: 0}``. +/// :rtype: dict +/// +/// .. [1] "Efficient Algorithms for Finding Maximum Matching in Graphs", +/// Zvi Galil, ACM Computing Surveys, 1986. +/// +#[pyfunction(max_cardinality = "false", default_weight = 1)] +#[text_signature = "(graph, /, max_cardinality=False, weight_fn=None, default_weight=1)"] +pub fn max_weight_matching( + py: Python, + graph: &graph::PyGraph, + max_cardinality: bool, + weight_fn: Option, + default_weight: i128, +) -> PyResult> { + max_weight_matching::max_weight_matching( + py, + graph, + max_cardinality, + weight_fn, + default_weight, + ) +} + /// Compute the strongly connected components for a directed graph /// /// This function is implemented using Kosaraju's algorithm @@ -2644,6 +2700,7 @@ fn retworkx(py: Python<'_>, m: &PyModule) -> PyResult<()> { m.add_wrapped(wrap_pyfunction!(digraph_find_cycle))?; m.add_wrapped(wrap_pyfunction!(digraph_k_shortest_path_lengths))?; m.add_wrapped(wrap_pyfunction!(graph_k_shortest_path_lengths))?; + m.add_wrapped(wrap_pyfunction!(max_weight_matching))?; m.add_class::()?; m.add_class::()?; m.add_class::()?; diff --git a/src/max_weight_matching.rs b/src/max_weight_matching.rs new file mode 100644 index 000000000..7ffc709d9 --- /dev/null +++ b/src/max_weight_matching.rs @@ -0,0 +1,1410 @@ +// Licensed under the Apache License, Version 2.0 (the "License"); you may +// not use this file except in compliance with the License. You may obtain +// a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, WITHOUT +// WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the +// License for the specific language governing permissions and limitations +// under the License. + +// Needed to pass shared state between functions +// closures don't work because of recurssion +#![allow(clippy::too_many_arguments)] +// Allow single character names to match naming convention from +// paper +#![allow(clippy::many_single_char_names)] + +use std::cmp::max; +use std::mem; + +use hashbrown::HashMap; + +use petgraph::graph::NodeIndex; +use petgraph::visit::{EdgeRef, IntoEdgeReferences}; + +use crate::graph::PyGraph; +use pyo3::exceptions::PyException; +use pyo3::prelude::*; + +fn weight_callable( + py: Python, + weight: PyObject, + weight_fn: &Option, + default: i128, +) -> PyResult { + match weight_fn { + Some(weight_fn) => { + let res = weight_fn.call1(py, (weight,))?; + res.extract(py) + } + None => Ok(default), + } +} + +// An inner container for a bloossom Vec +// if exists is false it hasn't been initialized yet +struct BlossomList { + exists: bool, + nodes: Vec, +} + +/// Return 2 * slack of edge k (does not work inside blossoms). +fn slack( + edge_index: usize, + dual_var: &[i128], + edges: &[(usize, usize, i128)], +) -> i128 { + let (source_index, target_index, weight) = edges[edge_index]; + dual_var[source_index] + dual_var[target_index] - 2 * weight +} + +/// Generate the leaf vertices of a blossom. +fn blossom_leaves( + blossom: usize, + num_nodes: usize, + blossom_children: &[BlossomList], +) -> PyResult> { + let mut out_vec: Vec = Vec::new(); + if blossom < num_nodes { + out_vec.push(blossom); + } else { + let child_blossom = &blossom_children[blossom]; + if !child_blossom.exists { + return Err(PyException::new_err("child blossom not initialized")); + } + for c in &child_blossom.nodes { + let child = *c; + if child < num_nodes { + out_vec.push(child); + } else { + for v in blossom_leaves(child, num_nodes, blossom_children)? { + out_vec.push(v); + } + } + } + } + Ok(out_vec) +} + +/// Assign label t to the top-level blossom containing vertex w +/// and record the fact that w was reached through the edge with +/// remote endpoint p. +fn assign_label( + w: usize, + t: usize, + p: Option, + num_nodes: usize, + in_blossoms: &[usize], + labels: &mut Vec>, + label_ends: &mut Vec>, + best_edge: &mut Vec>, + queue: &mut Vec, + blossom_children: &[BlossomList], + blossom_base: &[Option], + endpoints: &[usize], + mate: &[Option], +) -> PyResult<()> { + let b = in_blossoms[w]; + assert!(labels[w] == Some(0) && labels[b] == Some(0)); + labels[w] = Some(t); + labels[b] = Some(t); + label_ends[b] = p; + label_ends[w] = p; + label_ends[b] = p; + best_edge[w] = None; + best_edge[b] = None; + if t == 1 { + // b became an S-vertex/blossom; add it(s verticies) to the queue + queue.append(&mut blossom_leaves(b, num_nodes, &blossom_children)?); + } else if t == 2 { + // b became a T-vertex/blossom; assign label S to its mate. + // (If b is a non-trivial blossom, its base is the only vertex + // with an external mate.) + let blossom_index: usize = b; + let base: usize = blossom_base[blossom_index].unwrap(); + assert!(mate[base].is_some()); + assign_label( + endpoints[mate[base].unwrap()], + 1, + match mate[base] { + Some(p) => Some(p ^ 1), + None => None, + }, + num_nodes, + in_blossoms, + labels, + label_ends, + best_edge, + queue, + blossom_children, + blossom_base, + endpoints, + mate, + )?; + } + Ok(()) +} + +/// Trace back from vertices v and w to discover either a new blossom +/// or an augmenting path. Return the base vertex of the new blossom or None. +fn scan_blossom( + node_a: usize, + node_b: usize, + in_blossoms: &[usize], + blossom_base: &[Option], + endpoints: &[usize], + labels: &mut Vec>, + label_ends: &[Option], + mate: &[Option], +) -> Option { + let mut v: Option = Some(node_a); + let mut w: Option = Some(node_b); + // Trace back from v and w, placing breadcrumbs as we go + let mut path: Vec = Vec::new(); + let mut base: Option = None; + while v.is_some() || w.is_some() { + // Look for a breadcrumb in v's blossom or put a new breadcrump + let mut blossom = in_blossoms[v.unwrap()]; + if labels[blossom].is_none() || labels[blossom].unwrap() & 4 > 0 { + base = blossom_base[blossom]; + break; + } + assert!(labels[blossom] == Some(1)); + path.push(blossom); + labels[blossom] = Some(5); + // Trace one step bacl. + assert!(label_ends[blossom] == mate[blossom_base[blossom].unwrap()]); + if label_ends[blossom].is_none() { + // The base of blossom is single; stop tracing this path + v = None; + } else { + let tmp = endpoints[label_ends[blossom].unwrap()]; + blossom = in_blossoms[tmp]; + assert!(labels[blossom] == Some(2)); + // blossom is a T-blossom; trace one more step back. + assert!(label_ends[blossom].is_some()); + v = Some(endpoints[label_ends[blossom].unwrap()]); + } + // Swap v and w so that we alternate between both paths. + if w.is_some() { + mem::swap(&mut v, &mut w); + } + } + // Remvoe breadcrumbs. + for blossom in path { + labels[blossom] = Some(1); + } + // Return base vertex, if we found one. + base +} + +/// Construct a new blossom with given base, containing edge k which +/// connects a pair of S vertices. Label the new blossom as S; set its dual +/// variable to zero; relabel its T-vertices to S and add them to the queue. +fn add_blossom( + base: usize, + edge: usize, + blossom_children: &mut Vec, + num_nodes: usize, + edges: &[(usize, usize, i128)], + in_blossoms: &mut Vec, + dual_var: &mut Vec, + labels: &mut Vec>, + label_ends: &mut Vec>, + best_edge: &mut Vec>, + queue: &mut Vec, + blossom_base: &mut Vec>, + endpoints: &[usize], + blossom_endpoints: &mut Vec, + unused_blossoms: &mut Vec, + blossom_best_edges: &mut Vec, + blossom_parents: &mut Vec>, + neighbor_endpoints: &[Vec], + mate: &[Option], +) -> PyResult<()> { + let (mut v, mut w, _weight) = edges[edge]; + let blossom_b = in_blossoms[base]; + let mut blossom_v = in_blossoms[v]; + let mut blossom_w = in_blossoms[w]; + // Create blossom + let blossom = unused_blossoms.pop().unwrap(); + blossom_base[blossom] = Some(base); + blossom_parents[blossom] = None; + blossom_parents[blossom_b] = Some(blossom); + // Make list of sub-blossoms and their interconnecting edge endpoints. + blossom_children[blossom].exists = true; + let mut path: Vec = Vec::new(); + blossom_endpoints[blossom].exists = true; + let mut local_endpoints: Vec = Vec::new(); + // Trace back from blossom_v to base. + while blossom_v != blossom_b { + // Add blossom_v to the new blossom + blossom_parents[blossom_v] = Some(blossom); + path.push(blossom_v); + let blossom_v_endpoint_label = label_ends[blossom_v].unwrap(); + local_endpoints.push(blossom_v_endpoint_label); + assert!( + labels[blossom_v] == Some(2) + || (labels[blossom_v] == Some(1) + && label_ends[blossom_v] + == mate[blossom_base[blossom_v].unwrap()]) + ); + // Traceone step back. + assert!(label_ends[blossom_v].is_some()); + v = endpoints[blossom_v_endpoint_label]; + blossom_v = in_blossoms[v]; + } + // Reverse lists, add endpoint that connects the pair of S vertices. + path.push(blossom_b); + path.reverse(); + local_endpoints.reverse(); + local_endpoints.push(2 * edge); + // Trace back from w to base. + while blossom_w != blossom_b { + // Add blossom_w to the new blossom + blossom_parents[blossom_w] = Some(blossom); + path.push(blossom_w); + let blossom_w_endpoint_label = label_ends[blossom_w].unwrap(); + local_endpoints.push(blossom_w_endpoint_label ^ 1); + assert!( + labels[blossom_w] == Some(2) + || (labels[blossom_w] == Some(1) + && label_ends[blossom_w] + == mate[blossom_base[blossom_w].unwrap()]) + ); + // Trace one step back + assert!(label_ends[blossom_w].is_some()); + w = endpoints[blossom_w_endpoint_label]; + blossom_w = in_blossoms[w]; + } + // Set label to S. + assert!(labels[blossom_b] == Some(1)); + labels[blossom] = Some(1); + label_ends[blossom] = label_ends[blossom_b]; + // Set dual variable to 0 + dual_var[blossom] = 0; + // Relabel vertices + for node in blossom_leaves(blossom, num_nodes, &blossom_children)? { + if labels[in_blossoms[node]] == Some(2) { + // This T-vertex now turns into an S-vertex because it becomes + // part of an S-blossom; add it to the queue + queue.push(node); + } + in_blossoms[node] = blossom; + } + // Compute blossom_best_edges[blossom] + let mut best_edge_to: Vec> = vec![None; 2 * num_nodes]; + for bv in path { + // This sub-blossom does not have a list of least-slack edges; + // get the information from the vertices. + let nblists: Vec> = if !blossom_best_edges[bv].exists { + let mut tmp: Vec> = Vec::new(); + for node in blossom_leaves(bv, num_nodes, &blossom_children)? { + tmp.push( + neighbor_endpoints[node].iter().map(|p| p / 2).collect(), + ); + } + tmp + } else { + // Walk this sub-blossom's least-slack edges. + vec![blossom_best_edges[bv].nodes.clone()] + }; + for nblist in nblists { + for edge_index in nblist { + let (mut i, mut j, _edge_weight) = edges[edge_index]; + if in_blossoms[j] == blossom { + mem::swap(&mut i, &mut j); + } + let blossom_j = in_blossoms[j]; + if blossom_j != blossom + && labels[blossom_j] == Some(1) + && (best_edge_to[blossom_j].is_none() + || slack(edge_index, &dual_var, &edges) + < slack( + best_edge_to[blossom_j].unwrap(), + &dual_var, + &edges, + )) + { + best_edge_to[blossom_j] = Some(edge_index); + } + } + } + // Forget about least-slack edges of the sub-blossom. + blossom_best_edges[bv].exists = false; + blossom_best_edges[bv].nodes = Vec::new(); + best_edge[bv] = None; + } + blossom_best_edges[blossom].nodes = best_edge_to + .iter() + .filter(|x| x.is_some()) + .map(|x| x.unwrap()) + .collect(); + blossom_best_edges[blossom].exists = true; + //select best_edge[blossom] + best_edge[blossom] = None; + for edge_index in &blossom_best_edges[blossom].nodes { + if best_edge[blossom].is_none() + || slack(*edge_index, &dual_var, &edges) + < slack(best_edge[blossom].unwrap(), &dual_var, &edges) + { + best_edge[blossom] = Some(*edge_index); + } + } + Ok(()) +} + +/// Expand the given top level blossom +fn expand_blossom( + blossom: usize, + end_stage: bool, + num_nodes: usize, + blossom_children: &mut Vec, + in_blossoms: &mut Vec, + dual_var: &[i128], + labels: &mut Vec>, + label_ends: &mut Vec>, + best_edge: &mut Vec>, + queue: &mut Vec, + blossom_base: &[Option], + endpoints: &[usize], + mate: &[Option], + blossom_endpoints: &mut Vec, + allowed_edge: &mut Vec, + unused_blossoms: &mut Vec, +) -> PyResult<()> { + if !blossom_children[blossom].exists { + panic!("blossom children not defined"); + } + // convert sub-blossoms into top-level blossoms. + for s in blossom_children[blossom].nodes.clone() { + if s < num_nodes { + in_blossoms[s] = s + } else if end_stage && dual_var[s] == 0 { + // Recursively expand this sub-blossom + expand_blossom( + s, + end_stage, + num_nodes, + blossom_children, + in_blossoms, + dual_var, + labels, + label_ends, + best_edge, + queue, + blossom_base, + endpoints, + mate, + blossom_endpoints, + allowed_edge, + unused_blossoms, + )?; + } else { + for v in blossom_leaves(s, num_nodes, blossom_children)? { + in_blossoms[v] = s; + } + } + } + // if we expand a T-blossom during a stage, its a sub-blossoms must be + // relabeled + if !end_stage && labels[blossom] == Some(2) { + // start at the sub-blossom through which the expanding blossom + // obtained its label, and relabel sub-blossoms until we reach the + // base. + assert!(label_ends[blossom].is_some()); + let entry_child = + in_blossoms[endpoints[label_ends[blossom].unwrap() ^ 1]]; + // Decied in which direction we will go around the blossom. + let mut j = blossom_children[blossom] + .nodes + .iter() + .position(|x| *x == entry_child) + .unwrap(); + let j_step: i64; + let endpoint_trick: usize; + if j & 1 > 0 { + // Start index is odd; go forward and wrap + j -= blossom_children[blossom].nodes.len(); + j_step = 1; + endpoint_trick = 0; + } else { + // start index is even; go backward + j_step = -1; + endpoint_trick = 1; + } + // Move along the blossom until we get to the base + let mut p = label_ends[blossom].unwrap(); + while j != 0 { + // Relabel the T-sub-blossom. + labels[endpoints[p ^ 1]] = Some(0); + labels[endpoints[blossom_endpoints[blossom].nodes + [j - endpoint_trick] + ^ endpoint_trick + ^ 1]] = Some(0); + assign_label( + endpoints[p ^ 1], + 2, + Some(p), + num_nodes, + in_blossoms, + labels, + label_ends, + best_edge, + queue, + blossom_children, + blossom_base, + endpoints, + mate, + )?; + // Step to the next S-sub-blossom and note it's forwward endpoint. + allowed_edge + [blossom_endpoints[blossom].nodes[j - endpoint_trick] / 2] = + true; + if j_step > 0 { + j += j_step as usize; + } else { + j -= j_step.abs() as usize; + } + p = blossom_endpoints[blossom].nodes[j - endpoint_trick] + ^ endpoint_trick; + // Step to the next T-sub-blossom. + allowed_edge[p / 2] = true; + if j_step > 0 { + j += j_step as usize; + } else { + j -= j_step.abs() as usize; + } + } + // Relabel the base T-sub-blossom WITHOUT stepping through to + // its mate (so don't call assign_label()) + let blossom_v = blossom_children[blossom].nodes[j]; + labels[endpoints[p ^ 1]] = Some(2); + labels[blossom_v] = Some(2); + label_ends[endpoints[p ^ 1]] = Some(p); + label_ends[blossom_v] = Some(p); + best_edge[blossom_v] = None; + // Continue along the blossom until we get back to entry_child + while blossom_children[blossom].nodes[j] != entry_child { + // Examine the vertices of the sub-blossom to see whether + // it is reachable from a neighboring S-vertex outside the + // expanding blososm. + let bv = blossom_children[blossom].nodes[j]; + if labels[bv] == Some(1) { + // This sub-blossom just got label S through one of its + // neighbors; leave it + if j_step > 0 { + j += j_step as usize; + } else { + j -= j_step.abs() as usize; + } + continue; + } + let mut v: usize = 0; + for tmp in blossom_leaves(bv, num_nodes, blossom_children)? { + v = tmp; + if labels[v].unwrap() != 0 { + break; + } + } + // If the sub-blossom contains a reachable vertex, assign label T + // to the sub-blossom + if labels[v] != Some(0) { + assert!(labels[v] == Some(2)); + assert!(in_blossoms[v] == bv); + labels[v] = Some(0); + labels[endpoints[mate[blossom_base[bv].unwrap()].unwrap()]] = + Some(0); + assign_label( + v, + 2, + label_ends[v], + num_nodes, + in_blossoms, + labels, + label_ends, + best_edge, + queue, + blossom_children, + blossom_base, + endpoints, + mate, + )?; + } + if j_step > 0 { + j += j_step as usize; + } else { + j -= j_step.abs() as usize; + } + } + } + // Recycle the blossom number. + labels[blossom] = None; + label_ends[blossom] = None; + blossom_children[blossom].exists = false; + blossom_children[blossom].nodes.clear(); + blossom_endpoints[blossom].exists = false; + blossom_endpoints[blossom].nodes.clear(); + best_edge[blossom] = None; + unused_blossoms.push(blossom); + Ok(()) +} + +/// Swap matched/unmatched edges over an alternating path through blossom b +/// between vertex v and the base vertex. Keep blossom bookkeeping consistent. +fn augment_blossom( + blossom: usize, + node: usize, + num_nodes: usize, + blossom_parents: &[Option], + endpoints: &[usize], + blossom_children: &mut Vec, + blossom_endpoints: &mut Vec, + blossom_base: &mut Vec>, + mate: &mut Vec>, +) { + // Bubble up through the blossom tree from vertex v to an immediate + // sub-blossom of b. + let mut tmp = node; + while blossom_parents[tmp].unwrap() != blossom { + tmp = blossom_parents[tmp].unwrap(); + } + // Recursively deal with the first sub-blossom. + if tmp >= num_nodes { + augment_blossom( + tmp, + node, + num_nodes, + blossom_parents, + endpoints, + blossom_children, + blossom_endpoints, + blossom_base, + mate, + ); + } + // Decide in which direction we will go round the blossom. + let i = blossom_children[blossom] + .nodes + .iter() + .position(|x| *x == tmp) + .unwrap(); + let mut j = i; + let j_step: i64; + let endpoint_trick: usize; + if i & 1 > 0 { + // start index is odd; go forward and wrap. + j -= blossom_children[blossom].nodes.len(); + j_step = 1; + endpoint_trick = 0; + } else { + // Start index is even; go backward. + j_step = -1; + endpoint_trick = 1; + } + // Move along the blossom until we get to the base. + while j != 0 { + // Step to the next sub-blossom and augment it recursively. + if j_step > 0 { + j += j_step as usize; + } else { + j -= j_step.abs() as usize; + } + tmp = blossom_children[blossom].nodes[j]; + let p = blossom_endpoints[blossom].nodes[j - endpoint_trick] + ^ endpoint_trick; + if tmp > num_nodes { + augment_blossom( + tmp, + endpoints[p], + num_nodes, + blossom_parents, + endpoints, + blossom_children, + blossom_endpoints, + blossom_base, + mate, + ); + } + if j_step > 0 { + j += j_step as usize; + } else { + j -= j_step.abs() as usize; + } + tmp = blossom_children[blossom].nodes[j]; + if tmp > num_nodes { + augment_blossom( + tmp, + endpoints[p ^ 1], + num_nodes, + blossom_parents, + endpoints, + blossom_children, + blossom_endpoints, + blossom_base, + mate, + ); + } + // Match the edge connecting those sub-blossoms. + mate[endpoints[p]] = Some(p ^ 1); + mate[endpoints[p ^ 1]] = Some(p); + } + // Rotate the list of sub-blossoms to put the new base at the front. + let mut children: Vec = + blossom_children[blossom].nodes[i..].to_vec(); + children.extend_from_slice(&blossom_children[blossom].nodes[..i]); + blossom_children[blossom].exists = true; + blossom_children[blossom].nodes = children; + let mut endpoint: Vec = + blossom_endpoints[blossom].nodes[i..].to_vec(); + endpoint.extend_from_slice(&blossom_endpoints[blossom].nodes[..i]); + blossom_endpoints[blossom].exists = true; + blossom_endpoints[blossom].nodes = endpoint; + blossom_base[blossom] = blossom_base[blossom_children[blossom].nodes[0]]; + assert!(blossom_base[blossom] == Some(node)); +} + +/// Swap matched/unmatched edges over an alternating path between two +/// single vertices. The augmenting path runs through edge k, which +/// connects a pair of S vertices. +fn augment_matching( + edge: usize, + num_nodes: usize, + edges: &[(usize, usize, i128)], + in_blossoms: &[usize], + labels: &[Option], + label_ends: &[Option], + blossom_parents: &[Option], + endpoints: &[usize], + blossom_children: &mut Vec, + blossom_endpoints: &mut Vec, + blossom_base: &mut Vec>, + mate: &mut Vec>, +) { + let (v, w, _weight) = edges[edge]; + for (s_ref, p_ref) in [(v, 2 * edge + 1), (w, 2 * edge)].iter() { + // Match vertex s to remote endpoint p. Then trace back from s + // until we find a single vertex, swapping matched and unmatched + // edges as we go. + let mut s: usize = *s_ref; + let mut p: usize = *p_ref; + loop { + let blossom_s = in_blossoms[s]; + assert!(labels[blossom_s] == Some(1)); + assert!( + label_ends[blossom_s] == mate[blossom_base[blossom_s].unwrap()] + ); + // Augment through the S-blossom from s to base. + if blossom_s >= num_nodes { + augment_blossom( + blossom_s, + s, + num_nodes, + blossom_parents, + endpoints, + blossom_children, + blossom_endpoints, + blossom_base, + mate, + ); + } + // Update mate[s] + mate[s] = Some(p); + // Trace one step back. + if label_ends[blossom_s].is_none() { + // Reached a single vertex; stop + break; + } + let t = endpoints[label_ends[blossom_s].unwrap()]; + let blossom_t = in_blossoms[t]; + assert!(labels[blossom_t] == Some(2)); + // Trace one step back + assert!(label_ends[blossom_t].is_some()); + s = endpoints[label_ends[blossom_t].unwrap()]; + let j = endpoints[label_ends[blossom_t].unwrap() ^ 1]; + // Augment through the T-blossom from j to base. + assert!(blossom_base[blossom_t] == Some(t)); + if blossom_t >= num_nodes { + augment_blossom( + blossom_t, + j, + num_nodes, + blossom_parents, + endpoints, + blossom_children, + blossom_endpoints, + blossom_base, + mate, + ); + } + // Update mate[j] + mate[j] = label_ends[blossom_t]; + // Keep the opposite endpoint; + // it will be assigned to mate[s] in the next step. + p = label_ends[blossom_t].unwrap() ^ 1; + } + } +} + +///// Swap matched/unmatched edges over an alternating path between two +///// single vertices. The augmenting path runs through the edge, which +///// connects a pair of S vertices. +//fn verify_optimum( +// max_cardinality: bool, +// num_nodes: usize, +// num_edges: usize, +// edges: &[(usize, usize, i128)], +// endpoints: &[usize], +// dual_var: &[i128], +// blossom_parents: &[Option], +// blossom_endpoints: &[BlossomList], +// blossom_base: &[Option], +// mate: &[Option], +//) -> bool { +// let dual_var_node_min: i128 = *dual_var[..num_nodes].iter().min().unwrap(); +// let node_dual_offset: i128 = if max_cardinality { +// // Vertices may have negative dual; +// // find a constant non-negative number to add to all vertex duals. +// max(0, -dual_var_node_min) +// } else { +// 0 +// }; +// if dual_var_node_min + node_dual_offset < 0 { +// return false; +// } +// if *dual_var[num_nodes..].iter().min().unwrap() < 0 { +// return false; +// } +// // 0. all edges have non-negative slack and +// // 1. all matched edges have zero slack; +// for (edge, (i, j, weight)) in edges.iter().enumerate().take(num_edges) { +// let mut s = dual_var[*i] + dual_var[*j] - 2 * weight; +// let mut i_blossoms: Vec = vec![*i]; +// let mut j_blossoms: Vec = vec![*j]; +// while blossom_parents[*i_blossoms.last().unwrap()].is_some() { +// i_blossoms +// .push(blossom_parents[*i_blossoms.last().unwrap()].unwrap()); +// } +// while blossom_parents[*j_blossoms.last().unwrap()].is_some() { +// j_blossoms +// .push(blossom_parents[*j_blossoms.last().unwrap()].unwrap()); +// } +// i_blossoms.reverse(); +// j_blossoms.reverse(); +// for (blossom_i, blossom_j) in i_blossoms.iter().zip(j_blossoms.iter()) { +// if blossom_i != blossom_j { +// break; +// } +// s += 2 * dual_var[*blossom_i]; +// } +// if s < 0 { +// return false; +// } +// if mate[*i].unwrap() / 2 == edge || mate[*j].unwrap() / 2 == edge { +// if mate[*i].unwrap() / 2 != edge || mate[*j].unwrap() / 2 != edge { +// return false; +// } +// if s == 0 { +// return false; +// } +// } +// } +// // 2. all single vertices have zero dual value; +// for node in 0..num_nodes { +// if mate[node].is_none() && dual_var[node] + node_dual_offset != 0 { +// return false; +// } +// } +// // 3. all blossoms with positive dual value are full. +// for blossom in num_nodes..2 * num_nodes { +// if blossom_base[blossom].is_some() && dual_var[blossom] > 0 { +// if blossom_endpoints[blossom].nodes.len() % 2 != 1 { +// return false; +// } +// for p in [ +// blossom_endpoints[blossom].nodes[1], +// blossom_endpoints[blossom].nodes[2], +// ] +// .iter() +// { +// if mate[endpoints[*p]].unwrap() != *p ^ 1 { +// return false; +// } +// if mate[endpoints[*p ^ 1]].unwrap() != *p { +// return false; +// } +// } +// } +// } +// true +//} + +/// Compute a maximum-weighted matching in the general undirected weighted +/// graph given by "edges". If `max_cardinality` is true, only +/// maximum-cardinality matchings are considered as solutions. +/// +/// The algorithm is based on "Efficient Algorithms for Finding Maximum +/// Matching in Graphs" by Zvi Galil, ACM Computing Surveys, 1986. +/// +/// Based on networkx implementation +/// https://github.com/networkx/networkx/blob/3351206a3ce5b3a39bb2fc451e93ef545b96c95b/networkx/algorithms/matching.py +/// +/// With reference to the standalone protoype implementation from: +/// http://jorisvr.nl/article/maximum-matching +/// +/// http://jorisvr.nl/files/graphmatching/20130407/mwmatching.py +/// +/// The function takes time O(n**3) +pub fn max_weight_matching( + py: Python, + graph: &PyGraph, + max_cardinality: bool, + weight_fn: Option, + default_weight: i128, +) -> PyResult> { + let mut out_map: HashMap = HashMap::new(); + let num_edges = graph.graph.edge_count(); + let num_nodes = graph.graph.node_count(); + // Exit fast for graph without edges + if num_edges == 0 { + return Ok(out_map); + } + // Node indicies in the PyGraph may not be contiguous however the + // algorithm operates on contiguous indices 0..num_nodes node_map maps the + // algorithm index to the PyGraph's index + let node_map: HashMap = graph + .graph + .node_indices() + .enumerate() + .map(|(index, node_index)| (node_index, index)) + .collect(); + let mut edges: Vec<(usize, usize, i128)> = Vec::with_capacity(num_edges); + let mut max_weight: i128 = 0; + for edge in graph.graph.edge_references() { + let edge_weight: i128 = weight_callable( + py, + edge.weight().clone_ref(py), + &weight_fn, + default_weight, + )?; + if edge_weight > max_weight { + max_weight = edge_weight; + }; + edges.push(( + node_map[&edge.source()], + node_map[&edge.target()], + edge_weight, + )); + } + // If p is an edge endpoint + // endpoint[p] is the node index to which endpoint p is attached + let endpoints: Vec = (0..2 * num_edges) + .map(|endpoint| { + let edge_tuple = edges[endpoint / 2]; + let out_value: usize; + if endpoint % 2 == 0 { + out_value = edge_tuple.0; + } else { + out_value = edge_tuple.1; + } + out_value + }) + .collect(); + // If v is a node/vertex + // neighbor_endpoints[v] is the list of remote endpoints of the edges + // attached to v. + // Not modified by algorithm (only mut to initially construct contents). + let mut neighbor_endpoints: Vec> = + (0..num_nodes).map(|_| Vec::new()).collect(); + for edge in 0..num_edges { + neighbor_endpoints[edges[edge].0].push(2 * edge + 1); + neighbor_endpoints[edges[edge].1].push(2 * edge); + } + // If v is a vertex, + // mate[v] is the remote endpoint of its matched edge, or None if it is + // single (i.e. endpoint[mate[v]] is v's partner vertex). + // Initially all vertices are single; updated during augmentation. + let mut mate: Vec> = vec![None; num_nodes]; + + // If b is a top-level blossom + // label[b] is 0 if b is unlabeled (free); + // 1 if b is a S-vertex/blossom; + // 2 if b is a T-vertex/blossom; + // The label of a vertex/node is found by looking at the label of its + // top-level containing blossom + // If v is a node/vertex inside a T-Blossom, + // label[v] is LabelType::T_Vertex if and only if v is reachable from + // an S_Vertex outside the blossom. + let mut labels: Vec>; + // If b is a labeled top-level blossom, + // label_ends[b] is the remote endpoint of the edge through which b is + // obtained. Its label, or None if b's base vertex is single. + // If v is a vertex inside a T-blossom and label[v] == 2, + // label_ends[v] is the remote endpoint of the edge through which v is + // reachable from outside the blossom + let mut label_ends: Vec> = vec![None; 2 * num_nodes]; + // If v is a vertex/node + // in_blossoms[v] is the top-level blossom to which v belongs. + // If v is a top-level vertex, v is itself a blossom (a trivial blossom) + // and in_blossoms[v] == v. + // Initially all nodes are top-level trivial blossoms. + let mut in_blossoms: Vec = graph + .graph + .node_indices() + .map(|node| node.index()) + .collect(); + // if b is a sub-blossom + // blossom_parents[b] is its immediate parent + // If b is a top-level blossom, blossom_parents[b] is None + let mut blossom_parents: Vec> = vec![None; 2 * num_nodes]; + // If b is a non-trivial (sub-)blossom, + // blossom_children[b] is an ordered list of its sub-blossoms, starting with + // the base and going round the blossom. + let mut blossom_children: Vec = (0..2 * num_nodes) + .map(|_| BlossomList { + exists: false, + nodes: Vec::new(), + }) + .collect(); + // If b is a (sub-)blossom, + // blossombase[b] is its base VERTEX (i.e. recursive sub-blossom). + let mut blossom_base: Vec> = + (0..num_nodes).map(Some).collect(); + blossom_base.append(&mut vec![None; num_nodes]); + // If b is a non-trivial (sub-)blossom, + // blossomendps[b] is a list of endpoints on its connecting edges, + // such that blossomendps[b][i] is the local endpoint of blossomchilds[b][i] + // on the edge that connects it to blossomchilds[b][wrap(i+1)]. + let mut blossom_endpoints: Vec = (0..2 * num_nodes) + .map(|_| BlossomList { + exists: false, + nodes: Vec::new(), + }) + .collect(); + // If v is a free vertex (or an unreached vertex inside a T-blossom), + // best_edge[v] is the edge to an S-vertex with least slack, + // or -1 if there is no such edge. + // If b is a (possibly trivial) top-level S-blossom, + // best_edge[b] is the least-slack edge to a different S-blossom, + // or -1 if there is no such edge. + let mut best_edge: Vec>; + // If b is a non-trivial top-level S-blossom, + // blossom_best_edges[b] is a list of least-slack edges to neighboring + // S-blossoms, or None if no such list has been computed yet. + // This is used for efficient computation of delta3. + let mut blossom_best_edges: Vec = (0..2 * num_nodes) + .map(|_| BlossomList { + exists: false, + nodes: Vec::new(), + }) + .collect(); + let mut unused_blossoms: Vec = (num_nodes..2 * num_nodes).collect(); + // If v is a vertex, + // dual_var[v] = 2 * u(v) where u(v) is the v's variable in the dual + // optimization problem (multiplication by two ensures integer values + // throughout the algorithm if all edge weights are integers). + // If b is a non-trivial blossom, + // dual_var[b] = z(b) where z(b) is b's variable in the dual optimization + // problem. + let mut dual_var: Vec = vec![max_weight; num_nodes]; + dual_var.append(&mut vec![0; num_nodes]); + // If allowed_edge[k] is true, edge k has zero slack in the optimization + // problem; if allowed_edge[k] is false, the edge's slack may or may not + // be zero. + let mut allowed_edge: Vec; + // Queue of newly discovered S-vertices + let mut queue: Vec = Vec::with_capacity(num_nodes); + + // Main loop: continue until no further improvement is possible + for _ in 0..num_nodes { + // Each iteration of this loop is a "stage". + // A stage finds an augmenting path and uses that to improve + // the matching. + + // Removal labels from top-level blossoms/vertices + labels = vec![Some(0); 2 * num_nodes]; + // Forget all about least-slack edges. + best_edge = vec![None; 2 * num_nodes]; + blossom_best_edges.splice( + num_nodes.., + (0..num_nodes).map(|_| BlossomList { + exists: false, + nodes: Vec::new(), + }), + ); + // Loss of labeling means that we can not be sure that currently + // allowable edges remain allowable througout this stage. + allowed_edge = vec![false; num_edges]; + // Make queue empty + queue.clear(); + // Label single blossom/vertices with S and put them in queue. + for v in 0..num_nodes { + if mate[v].is_none() && labels[in_blossoms[v]] == Some(0) { + assign_label( + v, + 1, + None, + num_nodes, + &in_blossoms, + &mut labels, + &mut label_ends, + &mut best_edge, + &mut queue, + &blossom_children, + &blossom_base, + &endpoints, + &mate, + )?; + } + } + // Loop until we succeed in augmenting the matching. + let mut augmented = false; + loop { + // Each iteration of this loop is a "substage". + // A substage tries to find an augmenting path; + // if found, the path is used to improve the matching and + // the stage ends. If there is no augmenting path, the + // primal-dual method is used to pump some slack out of + // the dual variables. + + // Continue labeling until all vertices which are reachable + // through an alternating path have got a label. + while !queue.is_empty() && !augmented { + // Take an S vertex from the queue + let v = queue.pop().unwrap(); + assert!(labels[in_blossoms[v]] == Some(1)); + + // Scan it's neighbors + for p in &neighbor_endpoints[in_blossoms[v]] { + let k = *p / 2; + let mut kslack = 0; + let w = endpoints[*p]; + // w is a neighbor of v + if in_blossoms[v] == in_blossoms[w] { + // this edge is internal to a blossom; ignore it + continue; + } + if !allowed_edge[k] { + kslack = slack(k, &dual_var, &edges); + if kslack <= 0 { + // edge k has zero slack -> it is allowable + allowed_edge[k] = true; + } + } + if allowed_edge[k] { + if labels[in_blossoms[w]] == Some(0) { + // (C1) w is a free vertex; + // label w with T and label its mate with S (R12). + assign_label( + w, + 2, + Some(*p ^ 1), + num_nodes, + &in_blossoms, + &mut labels, + &mut label_ends, + &mut best_edge, + &mut queue, + &blossom_children, + &blossom_base, + &endpoints, + &mate, + )?; + } else if labels[in_blossoms[w]] == Some(1) { + // (C2) w is an S-vertex (not in the same blossom); + // follow back-links to discover either an + // augmenting path or a new blossom. + let base = scan_blossom( + v, + w, + &in_blossoms, + &blossom_base, + &endpoints, + &mut labels, + &label_ends, + &mate, + ); + if base.is_some() { + // Found a new blossom; add it to the blossom + // bookkeeping and turn it into an S-blossom. + add_blossom( + base.unwrap(), + k, + &mut blossom_children, + num_nodes, + &edges, + &mut in_blossoms, + &mut dual_var, + &mut labels, + &mut label_ends, + &mut best_edge, + &mut queue, + &mut blossom_base, + &endpoints, + &mut blossom_endpoints, + &mut unused_blossoms, + &mut blossom_best_edges, + &mut blossom_parents, + &neighbor_endpoints, + &mate, + )?; + } else { + // Found an augmenting path; augment the + // matching and end this stage. + augment_matching( + k, + num_nodes, + &edges, + &in_blossoms, + &labels, + &label_ends, + &blossom_parents, + &endpoints, + &mut blossom_children, + &mut blossom_endpoints, + &mut blossom_base, + &mut mate, + ); + augmented = true; + break; + } + } else if labels[w] == Some(0) { + // w is inside a T-blossom, but w itself has not + // yet been reached from outside the blossom; + // mark it as reached (we need this to relabel + // during T-blossom expansion). + assert!(labels[in_blossoms[w]] == Some(2)); + labels[w] = Some(2); + label_ends[w] = Some(*p ^ 1); + } + } else if labels[in_blossoms[w]] == Some(1) { + // keep track of the least-slack non-allowable edge to + // a different S-blossom + let blossom = in_blossoms[v]; + if best_edge[blossom].is_none() + || kslack + < slack( + best_edge[blossom].unwrap(), + &dual_var, + &edges, + ) + { + best_edge[blossom] = Some(k); + } + } else if labels[w] == Some(0) { + // w is a free vertex (or an unreached vertex inside + // a T-blossom) but we can not reach it yet; + // keep track of the least-slack edge that reaches w. + if best_edge[w].is_none() + || kslack + < slack( + best_edge[w].unwrap(), + &dual_var, + &edges, + ) + { + best_edge[w] = Some(k) + } + } + } + } + if augmented { + break; + } + // There is no augmenting path under these constraints; + // compute delta and reduce slack in the optimization problem. + // (Note that our vertex dual variables, edge slacks and delta's + // are pre-multiplied by two.) + let mut delta_type = -1; + let mut delta: Option = None; + let mut delta_edge: Option = None; + let mut delta_blossom: Option = None; + + // Compute delta1: the minumum value of any vertex dual. + if !max_cardinality { + delta_type = 1; + delta = Some(*dual_var[..num_nodes].iter().min().unwrap()); + } + + // Compute delta2: the minimum slack on any edge between + // an S-vertex and a free vertex. + for v in 0..num_nodes { + if labels[in_blossoms[v]] == Some(0) && best_edge[v].is_some() { + let d = slack(best_edge[v].unwrap(), &dual_var, &edges); + if delta_type == -1 || d < delta.unwrap() { + delta = Some(d); + delta_type = 2; + delta_edge = best_edge[v]; + } + } + } + + // Compute delta3: half the minimum slack on any edge between a + // pair of S-blossoms + for blossom in 0..2 * num_nodes { + if blossom_parents[blossom].is_none() + && labels[blossom] == Some(1) + && best_edge[blossom].is_some() + { + let kslack = + slack(best_edge[blossom].unwrap(), &dual_var, &edges); + assert!(kslack % 2 == 0); + let d = Some(kslack / 2); + if delta_type == -1 || d < delta { + delta = d; + delta_type = 3; + delta_edge = best_edge[blossom]; + } + } + } + + // Compute delta4: minimum z variable of any T-blossom + for blossom in num_nodes..2 * num_nodes { + if blossom_base[blossom].is_some() + && blossom_parents[blossom].is_none() + && labels[blossom] == Some(2) + && (delta_type == -1 || dual_var[blossom] < delta.unwrap()) + { + delta = Some(dual_var[blossom]); + delta_type = 4; + delta_blossom = Some(blossom); + } + } + if delta_type == -1 { + //No further improvement possible; max-cardinality optimum + // reached. Do a final delta update to make the optimum + // verifyable + assert!(max_cardinality); + delta_type = 1; + delta = + Some(max(0, *dual_var[..num_nodes].iter().min().unwrap())); + } + + // Update dual variables according to delta. + for v in 0..num_nodes { + if labels[in_blossoms[v]] == Some(1) { + // S-vertex: 2*u = 2*u - 2*delta + dual_var[v] -= delta.unwrap(); + } else if labels[in_blossoms[v]] == Some(2) { + // T-vertex: 2*u = 2*u + 2*delta + dual_var[v] += delta.unwrap(); + } + } + for b in num_nodes..2 * num_nodes { + if blossom_base[b].is_some() && blossom_parents[b].is_none() { + if labels[b] == Some(1) { + // top-level S-blossom: z = z + 2*delta + dual_var[b] += delta.unwrap(); + } else if labels[b] == Some(2) { + // top-level T-blossom: z = z - 2*delta + dual_var[b] -= delta.unwrap(); + } + } + } + // Take action at the point where minimum delta occured. + if delta_type == 1 { + // No further improvement possible; optimum reached + break; + } else if delta_type == 2 { + // Use the least-slack edge to continue the search. + allowed_edge[delta_edge.unwrap()] = true; + let (mut i, mut j, _weight) = edges[delta_edge.unwrap()]; + if labels[in_blossoms[i]] == Some(0) { + mem::swap(&mut i, &mut j); + } + assert!(labels[in_blossoms[i]] == Some(1)); + queue.push(i); + } else if delta_type == 3 { + // Use the least-slack edge to continue the search. + allowed_edge[delta_edge.unwrap()] = true; + let (i, _j, _weight) = edges[delta_edge.unwrap()]; + assert!(labels[in_blossoms[i]] == Some(1)); + queue.push(i); + } else if delta_type == 4 { + // Expand the least-z blossom + expand_blossom( + delta_blossom.unwrap(), + false, + num_nodes, + &mut blossom_children, + &mut in_blossoms, + &dual_var, + &mut labels, + &mut label_ends, + &mut best_edge, + &mut queue, + &blossom_base, + &endpoints, + &mate, + &mut blossom_endpoints, + &mut allowed_edge, + &mut unused_blossoms, + )?; + } + // end of this substage + } + // Stop when no more augment paths can be found + if !augmented { + break; + } + + // End of a stage; expand all S-blossoms which have a dual_var == 0 + for blossom in num_nodes..2 * num_nodes { + if blossom_parents[blossom].is_none() + && blossom_base[blossom].is_some() + && labels[blossom] == Some(1) + && dual_var[blossom] == 0 + { + expand_blossom( + blossom, + true, + num_nodes, + &mut blossom_children, + &mut in_blossoms, + &dual_var, + &mut labels, + &mut label_ends, + &mut best_edge, + &mut queue, + &blossom_base, + &endpoints, + &mate, + &mut blossom_endpoints, + &mut allowed_edge, + &mut unused_blossoms, + )?; + } + } + } + // if !verify_optimum( + // max_cardinality, + // num_nodes, + // num_edges, + // &edges, + // &endpoints, + // &dual_var, + // &blossom_parents, + // &blossom_endpoints, + // &blossom_base, + // &mate, + // ) { + // return Err(PyException::new_err("Solution found is not optimal")); + // } + + // Transform mate[] such that mate[v] is the vertex to which v is paired + // Also handle holes in node indices from graph removals by mapping + // linear index to node index. + let node_list: Vec = graph.graph.node_indices().collect(); + for (index, node) in mate.iter().enumerate() { + if node.is_some() { + out_map.insert( + node_list[index].index(), + node_list[endpoints[node.unwrap()]].index(), + ); + } + } + Ok(out_map) +} diff --git a/tests/test_max_weight_matching.py b/tests/test_max_weight_matching.py new file mode 100644 index 000000000..1759e6819 --- /dev/null +++ b/tests/test_max_weight_matching.py @@ -0,0 +1,303 @@ +# Licensed under the Apache License, Version 2.0 (the "License"); you may +# not use this file except in compliance with the License. You may obtain +# a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, WITHOUT +# WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the +# License for the specific language governing permissions and limitations +# under the License. + +# These tests are adapated from the networkx test cases: +# https://github.com/networkx/networkx/blob/3351206a3ce5b3a39bb2fc451e93ef545b96c95b/networkx/algorithms/tests/test_matching.py + +import unittest + +import retworkx + + +class TestMaxWeightMatching(unittest.TestCase): + + def test_empty_graph(self): + graph = retworkx.PyGraph() + self.assertEqual(retworkx.max_weight_matching(graph), {}) + + def test_single_edge(self): + graph = retworkx.PyGraph() + graph.add_nodes_from([0, 1]) + graph.add_edges_from([(0, 1, 1)]) + self.assertEqual(retworkx.max_weight_matching(graph), {0: 1, 1: 0}) + + def test_single_self_edge(self): + graph = retworkx.PyGraph() + graph.extend_from_weighted_edge_list([(0, 0, 100)]) + self.assertEqual(retworkx.max_weight_matching(graph), {}) + + def test_small_graph(self): + graph = retworkx.PyGraph() + graph.extend_from_weighted_edge_list([(1, 2, 10), (2, 3, 11)]) + self.assertEqual( + retworkx.max_weight_matching(graph, weight_fn=lambda x: x), + {3: 2, 2: 3}) + + def test_path_graph(self): + graph = retworkx.PyGraph() + graph.extend_from_weighted_edge_list( + [(1, 2, 5), (2, 3, 11), (3, 4, 5)]) + self.assertEqual( + retworkx.max_weight_matching(graph, weight_fn=lambda x: x), + {2: 3, 3: 2}) + self.assertEqual( + retworkx.max_weight_matching(graph, True, weight_fn=lambda x: x), + {1: 2, 2: 1, 3: 4, 4: 3}) + + def test_negative_weights(self): + graph = retworkx.PyGraph() + graph.extend_from_weighted_edge_list([ + (1, 2, 2), + (1, 3, -2), + (2, 3, 1), + (2, 4, -1), + (3, 4, -6), + ]) + self.assertEqual( + retworkx.max_weight_matching(graph, weight_fn=lambda x: x), + {1: 2, 2: 1}) + self.assertEqual( + retworkx.max_weight_matching(graph, True, weight_fn=lambda x: x), + {1: 3, 2: 4, 3: 1, 4: 2}) + + def test_s_blossom(self): + graph = retworkx.PyGraph() + graph.extend_from_weighted_edge_list([ + (0, 1, 8), + (0, 2, 9), + (1, 2, 10), + (2, 3, 7), + ]) + self.assertEqual( + retworkx.max_weight_matching(graph, weight_fn=lambda x: x), + {0: 1, 1: 0, 2: 3, 3: 2}) + graph.extend_weighted_edge_list([(0, 5, 5), (3, 4, 6)]) + self.assertEqual( + retworkx.max_weight_matching(graph, weight_fn=lambda x: x), + {0: 5, 1: 2, 2: 1, 3: 4, 4: 3, 5: 0}) + + def test_s_t_blossom(self): + graph = retworkx.PyGraph() + graph.extend_from_weighted_edge_list([ + (1, 2, 9), + (1, 3, 8), + (2, 3, 10), + (1, 4, 5), + (4, 5, 4), + (1, 6, 3), + ]) + self.assertEqual( + retworkx.max_weight_matching(graph, weight_fn=lambda x: x), + {1: 6, 2: 3, 3: 2, 4: 5, 5: 4, 6: 1}) + graph.remove_edge(1, 6) + graph.remove_edge(4, 5) + graph.extend_from_weighted_edge_list([(4, 5, 3), (1, 6, 4)]) + self.assertEqual( + retworkx.max_weight_matching(graph, weight_fn=lambda x: x), + {1: 6, 2: 3, 3: 2, 4: 5, 5: 4, 6: 1}) + graph.remove_edge(1, 6) + graph.add_edge(3, 6, 4) + self.assertEqual( + retworkx.max_weight_matching(graph, weight_fn=lambda x: x), + {1: 2, 2: 1, 3: 6, 4: 5, 5: 4, 6: 3}) + + def test_nested_s_blossom(self): + graph = retworkx.PyGraph() + graph.extend_from_weighted_edge_list([ + (1, 2, 9), + (1, 3, 9), + (2, 3, 10), + (2, 4, 8), + (3, 5, 8), + (4, 5, 10), + (5, 6, 6), + ]) + expected = {1: 3, 2: 4, 3: 1, 4: 2, 5: 6, 6: 5} + self.assertEqual( + retworkx.max_weight_matching(graph, weight_fn=lambda x: x), + expected) + + def test_nested_s_blossom_relabel(self): + graph = retworkx.PyGraph() + graph.extend_from_weighted_edge_list([ + (1, 2, 10), + (1, 7, 10), + (2, 3, 12), + (3, 4, 20), + (3, 5, 20), + (4, 5, 25), + (5, 6, 10), + (6, 7, 10), + (7, 8, 8), + ]) + self.assertEqual( + retworkx.max_weight_matching(graph, weight_fn=lambda x: x), + {1: 2, 2: 1, 3: 4, 4: 3, 5: 6, 6: 5, 7: 8, 8: 7}) + + def test_nested_s_blossom_expand(self): + graph = retworkx.PyGraph() + graph.extend_from_weighted_edge_list([ + (1, 2, 8), + (1, 3, 8), + (2, 3, 10), + (2, 4, 12), + (3, 5, 12), + (4, 5, 14), + (4, 6, 12), + (5, 7, 12), + (6, 7, 14), + (7, 8, 12), + ]) + self.assertEqual( + retworkx.max_weight_matching(graph, weight_fn=lambda x: x), + {1: 2, 2: 1, 3: 5, 4: 6, 5: 3, 6: 4, 7: 8, 8: 7}) + + def test_s_blossom_relabel_expand(self): + graph = retworkx.PyGraph() + graph.extend_from_weighted_edge_list([ + (1, 2, 23), + (1, 5, 22), + (1, 6, 15), + (2, 3, 25), + (3, 4, 22), + (4, 5, 25), + (4, 8, 14), + (5, 7, 13), + ]) + self.assertEqual( + retworkx.max_weight_matching(graph, weight_fn=lambda x: x), + {1: 6, 2: 3, 3: 2, 4: 8, 5: 7, 6: 1, 7: 5, 8: 4}) + + def test_nested_s_blossom_relabel_expand(self): + graph = retworkx.PyGraph() + graph.extend_from_weighted_edge_list([ + (1, 2, 19), + (1, 3, 20), + (1, 8, 8), + (2, 3, 25), + (2, 4, 18), + (3, 5, 18), + (4, 5, 13), + (4, 7, 7), + (5, 6, 7), + ]) + self.assertEqual( + retworkx.max_weight_matching(graph, weight_fn=lambda x: x), + {1: 8, 2: 3, 3: 2, 4: 7, 5: 6, 6: 5, 7: 4, 8: 1}) + + def test_blossom_relabel_multiple_paths(self): + graph = retworkx.PyGraph() + graph.extend_from_weighted_edge_list([ + (1, 2, 45), + (1, 5, 45), + (2, 3, 50), + (3, 4, 45), + (4, 5, 50), + (1, 6, 30), + (3, 9, 35), + (4, 8, 35), + (5, 7, 26), + (9, 10, 5), + ]) + self.assertEqual( + retworkx.max_weight_matching(graph, weight_fn=lambda x: x), + {1: 6, 2: 3, 3: 2, 4: 8, 5: 7, 6: 1, 7: 5, 8: 4, 9: 10, 10: 9}) + + def test_blossom_relabel_multiple_path_alternate(self): + graph = retworkx.PyGraph() + graph.extend_from_weighted_edge_list([ + (1, 2, 45), + (1, 5, 45), + (2, 3, 50), + (3, 4, 45), + (4, 5, 50), + (1, 6, 30), + (3, 9, 35), + (4, 8, 26), + (5, 7, 40), + (9, 10, 5), + ]) + self.assertEqual( + retworkx.max_weight_matching(graph, weight_fn=lambda x: x), + {1: 6, 2: 3, 3: 2, 4: 8, 5: 7, 6: 1, 7: 5, 8: 4, 9: 10, 10: 9}) + + def test_blossom_relabel_multiple_paths_least_slack(self): + graph = retworkx.PyGraph() + graph.extend_from_weighted_edge_list([ + (1, 2, 45), + (1, 5, 45), + (2, 3, 50), + (3, 4, 45), + (4, 5, 50), + (1, 6, 30), + (3, 9, 35), + (4, 8, 28), + (5, 7, 26), + (9, 10, 5), + ]) + self.assertEqual( + retworkx.max_weight_matching(graph, weight_fn=lambda x: x), + {1: 6, 2: 3, 3: 2, 4: 8, 5: 7, 6: 1, 7: 5, 8: 4, 9: 10, 10: 9}) + + def test_nested_blossom_expand_recursively(self): + graph = retworkx.PyGraph() + graph.extend_from_weighted_edge_list([ + (1, 2, 40), + (1, 3, 40), + (2, 3, 60), + (2, 4, 55), + (3, 5, 55), + (4, 5, 50), + (1, 8, 15), + (5, 7, 30), + (7, 6, 10), + (8, 10, 10), + (4, 9, 30), + ]) + self.assertEqual( + retworkx.max_weight_matching(graph, weight_fn=lambda x: x), + {1: 2, 2: 1, 3: 5, 4: 9, 5: 3, 6: 7, 7: 6, 8: 10, 9: 4, 10: 8}) + + def test_nested_blossom_augmented(self): + graph = retworkx.PyGraph() + graph.extend_from_weighted_edge_list([ + (1, 2, 45), + (1, 7, 45), + (2, 3, 50), + (3, 4, 45), + (4, 5, 95), + (4, 6, 94), + (5, 6, 94), + (6, 7, 50), + (1, 8, 30), + (3, 11, 35), + (5, 9, 36), + (7, 10, 26), + (11, 12, 5), + ]) + expected = { + 1: 8, + 2: 3, + 3: 2, + 4: 6, + 5: 9, + 6: 4, + 7: 10, + 8: 1, + 9: 5, + 10: 7, + 11: 12, + 12: 11, + } + self.assertEqual( + retworkx.max_weight_matching(graph, weight_fn=lambda x: x), + expected)