From 4f62a15c06f6d85aba30e82f07cd7f95eecf8c72 Mon Sep 17 00:00:00 2001 From: Joris van Rantwijk Date: Mon, 6 Feb 2023 17:39:03 +0100 Subject: [PATCH] Use (i, j) instead of (v, w) --- python/max_weight_matching.py | 589 +++++++++++++++++----------------- 1 file changed, 298 insertions(+), 291 deletions(-) diff --git a/python/max_weight_matching.py b/python/max_weight_matching.py index 31d6b28..93ba72e 100644 --- a/python/max_weight_matching.py +++ b/python/max_weight_matching.py @@ -34,8 +34,8 @@ def maximum_weight_matching( This function uses O(n + m) memory, where "m" is the number of edges. Parameters: - edges: List of edges, each edge specified as a tuple "(i, j, wt)" - where "i" and "j" are vertex indices and "wt" is the edge weight. + edges: List of edges, each edge specified as a tuple "(i, j, w)" + where "i" and "j" are vertex indices and "w" is the edge weight. Returns: List of pairs of matched vertex indices. @@ -79,7 +79,7 @@ def maximum_weight_matching( # Extract the final solution. pairs: list[tuple[int, int]] = [ - (i, j) for (i, j, _wt) in edges if matching.vertex_mate[i] == j] + (i, j) for (i, j, _w) in edges if matching.vertex_mate[i] == j] # Verify that the matching is optimal. # This only works reliably for integer weights. @@ -131,8 +131,8 @@ def adjust_weights_for_maximum_cardinality_matching( This function takes time O(m), where "m" is the number of edges. Parameters: - edges: List of edges, each edge specified as a tuple "(i, j, wt)" - where "i" and "j" are vertex indices and "wt" is the edge weight. + edges: List of edges, each edge specified as a tuple "(i, j, w)" + where "i" and "j" are vertex indices and "w" is the edge weight. Returns: List of edges with adjusted weights. If no adjustments are necessary, @@ -149,10 +149,10 @@ def adjust_weights_for_maximum_cardinality_matching( if not edges: return edges - num_vertex = 1 + max(max(i, j) for (i, j, _wt) in edges) + num_vertex = 1 + max(max(i, j) for (i, j, _w) in edges) - min_weight = min(wt for (_i, _j, wt) in edges) - max_weight = max(wt for (_i, _j, wt) in edges) + min_weight = min(w for (_i, _j, w) in edges) + max_weight = max(w for (_i, _j, w) in edges) weight_range = max_weight - min_weight # Do nothing if the weights already ensure a maximum-cardinality matching. @@ -171,7 +171,7 @@ def adjust_weights_for_maximum_cardinality_matching( assert delta >= 0 # Increase all edge weights by "delta". - return [(i, j, wt + delta) for (i, j, wt) in edges] + return [(i, j, w + delta) for (i, j, w) in edges] def _check_input_types(edges: list[tuple[int, int, int|float]]) -> None: @@ -181,8 +181,8 @@ def _check_input_types(edges: list[tuple[int, int, int|float]]) -> None: This function takes time O(m). Parameters: - edges: List of edges, each edge specified as a tuple "(i, j, wt)" - where "i" and "j" are edge indices and "wt" is the edge weight. + edges: List of edges, each edge specified as a tuple "(i, j, w)" + where "i" and "j" are edge indices and "w" is the edge weight. Raises: ValueError: If the input does not satisfy the constraints. @@ -198,7 +198,7 @@ def _check_input_types(edges: list[tuple[int, int, int|float]]) -> None: if (not isinstance(e, tuple)) or (len(e) != 3): raise TypeError("Each edge must be specified as a 3-tuple") - (i, j, wt) = e + (i, j, w) = e if (not isinstance(i, int)) or (not isinstance(j, int)): raise TypeError("Edge endpoints must be integers") @@ -206,17 +206,17 @@ def _check_input_types(edges: list[tuple[int, int, int|float]]) -> None: if (i < 0) or (j < 0): raise ValueError("Edge endpoints must be non-negative integers") - if not isinstance(wt, (int, float)): + if not isinstance(w, (int, float)): raise TypeError( "Edge weights must be integers or floating point numbers") - if isinstance(wt, float): - if not math.isfinite(wt): + if isinstance(w, float): + if not math.isfinite(w): raise ValueError("Edge weights must be finite numbers") # Check that this edge weight will not cause our dual variable # calculations to exceed the valid floating point range. - if wt > float_limit: + if w > float_limit: raise ValueError("Floating point edge weights must be" f" less than {float_limit:g}") @@ -228,15 +228,15 @@ def _check_input_graph(edges: list[tuple[int, int, int|float]]) -> None: This function takes time O(m * log(m)). Parameters: - edges: List of edges, each edge specified as a tuple "(i, j, wt)" - where "i" and "j" are edge indices and "wt" is the edge weight. + edges: List of edges, each edge specified as a tuple "(i, j, w)" + where "i" and "j" are edge indices and "w" is the edge weight. Raises: ValueError: If the input does not satisfy the constraints. """ # Check that the graph has no self-edges. - for (i, j, _wt) in edges: + for (i, j, _w) in edges: if i == j: raise ValueError("Self-edges are not supported") @@ -244,8 +244,7 @@ def _check_input_graph(edges: list[tuple[int, int, int|float]]) -> None: # Using a set() would be more straightforward, but the runtime bounds # of the Python set type are not clearly specified. # Sorting provides guaranteed O(m * log(m)) run time. - edge_endpoints = [((i, j) if (i < j) else (j, i)) - for (i, j, _wt) in edges] + edge_endpoints = [((i, j) if (i < j) else (j, i)) for (i, j, _w) in edges] edge_endpoints.sort() for i in range(len(edge_endpoints) - 1): @@ -271,17 +270,17 @@ class _GraphInfo: # Each edge is incident on two vertices. # Each edge also has a weight. # - # "edges[e] = (i, j, wt)" where + # "edges[e] = (i, j, w)" where # "e" is an edge index; # "i" and "j" are vertex indices of the incident vertices; - # "wt" is the edge weight. + # "w" is the edge weight. # # These data remain unchanged while the algorithm runs. self.edges: list[tuple[int, int, int|float]] = edges # num_vertex = the number of vertices. if edges: - self.num_vertex = 1 + max(max(i, j) for (i, j, _wt) in edges) + self.num_vertex = 1 + max(max(i, j) for (i, j, _w) in edges) else: self.num_vertex = 0 @@ -293,14 +292,14 @@ class _GraphInfo: # These data remain unchanged while the algorithm runs. self.adjacent_edges: list[list[int]] = [ [] for v in range(self.num_vertex)] - for (e, (i, j, _wt)) in enumerate(edges): + for (e, (i, j, _w)) in enumerate(edges): self.adjacent_edges[i].append(e) self.adjacent_edges[j].append(e) # Determine whether _all_ weights are integers. # In this case we can avoid floating point computations entirely. - self.integer_weights: bool = all(isinstance(wt, int) - for (_i, _j, wt) in edges) + self.integer_weights: bool = all(isinstance(w, int) + for (_i, _j, w) in edges) # TODO: @@ -316,17 +315,19 @@ class _GraphInfo: # - Is there a way to reduce allocations of tuples ? # (First profile if it saves any time, otherwise don't even bother.) # -# - Revert to using variable "w" for edge weight. -# - Avoid using "v" and "w" for vertex indices; instead use "i", "j", "k", or maybe "p", "q", "r". -# - Consistently use "b", "bi", "bj", "bk" for blossom index. -# And use "blossom", "iblossom", "jblossom" for the full blossom object. -# # - Maybe nice to abstract away the management of least-slack edges to separate classes. # # NOTE: It is possible for a blossom to have an unmatched base vertex. # Top-level blossoms with unmatched based vertex are marked S-blossom at the start of a stage. # Such blossoms can appear as the endpoints of an augmenting path. # +# NOTE: We use least-slack edge tracking between S-vertex and T-vertex +# to re-discover edges that can be used to relabel an unlabeled +# blossom after T-blossom expansion. +# This is a pretty good trick. Galil does it too. I probably did +# not understand this when I wrote the old code and figured out +# a different approach. +# # Proof that S-to-S edges have even slack when working with integer weights: # Edge slack is even iff indicent vertex duals are both odd or both even. # Unmatched vertices are always S vertices, therefore either all unmatched vertices are odd or all unmatched vertices are even. @@ -405,11 +406,11 @@ class _Blossom: self.subblossoms: list[int] = subblossoms # "edges" is a list of edges linking the sub-blossoms. - # Each edge is represented as an ordered pair "(v, w)" where "v" - # and "w" are vertex indices. + # Each edge is represented as an ordered pair "(i, j)" where "i" + # and "j" are vertex indices. # - # "edges[0] = (v, w)" where vertex "v" in "subblossoms[0]" is - # adjacent to vertex "w" in "subblossoms[1]", etc. + # "edges[0] = (i, j)" where vertex "i" in "subblossoms[0]" is + # adjacent to vertex "j" in "subblossoms[1]", etc. self.edges: list[tuple[int, int]] = edges # "base_vertex" is the vertex index of the base of the blossom. @@ -439,9 +440,9 @@ class _PartialMatching: # Each vertex is either single (unmatched) or matched to # another vertex. # - # If vertex "v" is matched to vertex "w", - # "vertex_mate[v] == w" and "vertex_mate[w] == v". - # If vertex "v" is unmatched, "vertex_mate[v] == -1". + # If vertex "i" is matched to vertex "j", + # "vertex_mate[i] == j" and "vertex_mate[j] == i". + # If vertex "i" is unmatched, "vertex_mate[i] == -1". # # Initially all vertices are unmatched. self.vertex_mate: list[int] = graph.num_vertex * [-1] @@ -465,22 +466,22 @@ class _PartialMatching: # List of currently unused blossom indices. self.unused_blossoms: list[int] = list( - range(graph.num_vertex, 2 * graph.num_vertex)) + reversed(range(graph.num_vertex, 2 * graph.num_vertex))) # TODO : we may need a list of top-level blossom indices # Every vertex is part of exactly one top-level blossom, # possibly a trivial blossom consisting of just that vertex. # - # "vertex_blossom[v]" is the index of the top-level blossom that - # contains vertex "v". - # "vertex_blossom[v] == v" if the "v" is a trivial top-level blossom. + # "vertex_blossom[i]" is the index of the top-level blossom that + # contains vertex "i". + # "vertex_blossom[i] == i" if the "i" is a trivial top-level blossom. # # Initially all vertices are top-level trivial blossoms. self.vertex_blossom: list[int] = list(range(graph.num_vertex)) # "blossom_parent[b]" is the index of the smallest blossom that - # contains blossom "b", or + # is a strict superset of blossom "b", or # "blossom_parent[b] == -1" if blossom "b" is a top-level blossom. # # Initially all vertices are trivial top-level blossoms. @@ -488,12 +489,12 @@ class _PartialMatching: # Every vertex has a variable in the dual LPP. # - # "dual_var_2x[v]" is 2 times the dual variable of "v". + # "dual_var_2x[i]" is 2 times the dual variable of vertex "i". # Multiplication by 2 ensures that the values are integers # if all edge weights are integers. # # Vertex duals are initialized to half the maximum edge weight. - max_weight = max(wt for (_i, _j, wt) in graph.edges) + max_weight = max(w for (_i, _j, w) in graph.edges) self.dual_var_2x: list[int|float] = graph.num_vertex * [max_weight] def edge_slack_2x(self, e: int) -> int|float: @@ -505,9 +506,9 @@ class _PartialMatching: Multiplication by 2 ensures that the return value is an integer if all edge weights are integers. """ - (i, j, wt) = self.graph.edges[e] + (i, j, w) = self.graph.edges[e] assert self.vertex_blossom[i] != self.vertex_blossom[j] - return self.dual_var_2x[i] + self.dual_var_2x[j] - 2 * wt + return self.dual_var_2x[i] + self.dual_var_2x[j] - 2 * w def get_blossom(self, b: int) -> _Blossom: """Return the Blossom instance for blossom index "b".""" @@ -560,17 +561,17 @@ class _StageData: # For each labeled blossom, we keep track of the edge that attaches # it to its alternating tree. # - # "blossom_link[b] = (v, w)" denotes the edge through which - # blossom "b" is attached to the alternating tree, where "v" and "w" - # are vertex indices and vertex "w" is contained in blossom "b". + # "blossom_link[b] = (i, j)" denotes the edge through which + # blossom "b" is attached to the alternating tree, where "i" and "j" + # are vertex indices and vertex "j" is contained in blossom "b". # # "blossom_link[b] = None" if "b" is the root of an alternating tree, # or if "b" is not a labeled, top-level blossom. self.blossom_link: list[Optional[tuple[int, int]]] = [ None for b in range(2 * num_vertex)] - # "vertex_best_edge[v]" is the edge index of the least-slack edge - # between "v" and any S-vertex, or -1 if no such edge has been found. + # "vertex_best_edge[i]" is the edge index of the least-slack edge + # between "i" and any S-vertex, or -1 if no such edge has been found. self.vertex_best_edge: list[int] = num_vertex * [-1] # For non-trivial top-level S-blossom "b", @@ -600,10 +601,10 @@ class _AlternatingPath(NamedTuple): def _trace_alternating_paths( matching: _PartialMatching, stage_data: _StageData, - v: int, - w: int + i: int, + j: int ) -> _AlternatingPath: - """Trace back through the alternating trees from vertices "v" and "w". + """Trace back through the alternating trees from vertices "i" and "j". If both vertices are part of the same alternating tree, this function discovers a new blossom. In this case it returns an alternating path @@ -626,44 +627,44 @@ def _trace_alternating_paths( blossom_marker = stage_data.blossom_marker marked_blossoms: list[int] = [] - # "vedges" is a list of edges used while tracing from "v". - # "wedges" is a list of edges used while tracing from "w". - vedges: list[tuple[int, int]] = [] - wedges: list[tuple[int, int]] = [] + # "iedges" is a list of edges used while tracing from "i". + # "jedges" is a list of edges used while tracing from "j". + iedges: list[tuple[int, int]] = [] + jedges: list[tuple[int, int]] = [] - # Pre-load the edge between "v" and "w" so it will end up in the right + # Pre-load the edge between "i" and "j" so it will end up in the right # place in the final path. - vedges.append((v, w)) + iedges.append((i, j)) - # Alternate between tracing the path from "v" and the path from "w". + # Alternate between tracing the path from "i" and the path from "j". # This ensures that the search time is bounded by the size of the # newly found blossom. first_common = -1 - while v != -1 or w != -1: + while i != -1 or j != -1: # Check if we found a common ancestor. - vblossom = vertex_blossom[v] - if blossom_marker[vblossom]: - first_common = vblossom + bi = vertex_blossom[i] + if blossom_marker[bi]: + first_common = bi break # Mark blossom as a potential common ancestor. - blossom_marker[vblossom] = True - marked_blossoms.append(vblossom) + blossom_marker[bi] = True + marked_blossoms.append(bi) # Track back through the link in the alternating tree. - link = blossom_link[vblossom] + link = blossom_link[bi] if link is None: # Reached the root of this alternating tree. - v = -1 + i = -1 else: - vedges.append(link) - v = link[0] + iedges.append(link) + i = link[0] - # Swap "v" and "w" to alternate between paths. - if w != -1: - v, w = w, v - vedges, wedges = wedges, vedges + # Swap "i" and "j" to alternate between paths. + if j != -1: + i, j = j, i + iedges, jedges = jedges, iedges # Remove all markers we placed. for b in marked_blossoms: @@ -671,13 +672,14 @@ def _trace_alternating_paths( # If we found a common ancestor, trim the paths so they end there. if first_common != -1: - assert vertex_blossom[vedges[-1][0]] == first_common - while wedges and (vertex_blossom[wedges[-1][0]] != first_common): - wedges.pop() + assert vertex_blossom[iedges[-1][0]] == first_common + while jedges and (vertex_blossom[jedges[-1][0]] != first_common): + jedges.pop() # Fuse the two paths. - path_edges = vedges[::-1] - path_edges.extend((q, p) for (p, q) in wedges) + # Flip the order of one path and the edge tuples in the other path + # to obtain a single continuous path with correctly ordered edge tuples. + path_edges = iedges[::-1] + [(j, i) for (i, j) in jedges] # Any S-to-S alternating path must have odd length. assert len(path_edges) % 2 == 1 @@ -705,12 +707,12 @@ def _make_blossom( assert len(path.edges) >= 3 # Construct the list of sub-blossoms (current top-level blossoms). - subblossoms = [matching.vertex_blossom[v] for (v, w) in path.edges] + subblossoms = [matching.vertex_blossom[i] for (i, j) in path.edges] # Check that the path is cyclic. # Note the path may not start and end with the same _vertex_, # but it must start and end in the same _blossom_. - subblossoms_next = [matching.vertex_blossom[w] for (v, w) in path.edges] + subblossoms_next = [matching.vertex_blossom[j] for (i, j) in path.edges] assert subblossoms[0] == subblossoms_next[-1] assert subblossoms[1:] == subblossoms_next[:-1] @@ -738,8 +740,8 @@ def _make_blossom( # to O(n**2) total time per stage. # This could be improved through a union-find datastructure, or # by somehow re-using the blossom index of the largest sub-blossom. - for v in matching.blossom_vertices(b): - matching.vertex_blossom[v] = b + for i in matching.blossom_vertices(b): + matching.vertex_blossom[i] = b # Assign label S to the new blossom. assert stage_data.blossom_label[base_blossom] == _LABEL_S @@ -772,8 +774,8 @@ def _make_blossom( # For sparse graphs, this could be improved by tracking # least-slack edges in a priority queue. best_edge_to_blossom: list[int] = (2 * num_vertex) * [-1] - best_slack_to_blossom: list[int|float] = [ - 0 for b in range(2 * num_vertex)] + zero_slack: int|float = 0 + best_slack_to_blossom: list[int|float] = (2 * num_vertex) * [zero_slack] # Add the least-slack edges of every S-sub-blossom. for sub in subblossoms: @@ -794,26 +796,28 @@ def _make_blossom( # Add edges to the temporary array. for e in sub_edge_set: - (i, j, _wt) = matching.graph.edges[e] - iblossom = matching.vertex_blossom[i] - jblossom = matching.vertex_blossom[j] - assert (iblossom == b) or (jblossom == b) - vblossom = jblossom if iblossom == b else iblossom + (i, j, _w) = matching.graph.edges[e] + bi = matching.vertex_blossom[i] + bj = matching.vertex_blossom[j] + assert (bi == b) or (bj == b) # Reject internal edges in this blossom. - if vblossom == b: + if bi == bj: continue + # Set bi = other blossom which is reachable through this edge. + bi = bj if (bi == b) else bi + # Reject edges that don't link to an S-blossom. - if stage_data.blossom_label[vblossom] != _LABEL_S: + if stage_data.blossom_label[bi] != _LABEL_S: continue # Keep only the least-slack edge to "vblossom". slack = matching.edge_slack_2x(e) - if ((best_edge_to_blossom[vblossom] == -1) - or (slack < best_slack_to_blossom[vblossom])): - best_edge_to_blossom[vblossom] = e - best_slack_to_blossom[vblossom] = slack + if ((best_edge_to_blossom[bi] == -1) + or (slack < best_slack_to_blossom[bi])): + best_edge_to_blossom[bi] = e + best_slack_to_blossom[bi] = slack # Extract a compact list of least-slack edge indices. # We can not keep the temporary array because that would blow up @@ -835,92 +839,93 @@ def _make_blossom( def _substage_assign_label_s( matching: _PartialMatching, stage_data: _StageData, - v: int, + i: int, ) -> None: - """Assign label S to the unlabeled blossom that contains vertex "v". + """Assign label S to the unlabeled blossom that contains vertex "i". - If vertex "v" is matched, it is attached to the alternating tree via its - matched edge. If vertex "v" is unmatched, it becomes the root of + If vertex "i" is matched, it is attached to the alternating tree via its + matched edge. If vertex "i" is unmatched, it becomes the root of an alternating tree. All vertices in the newly labeled S-blossom are added to the scan queue. Precondition: - - "v" is an unlabeled vertex, either unmatched - or matched to a T-vertex via a tight edge. + "i" is an unlabeled vertex, either unmatched or matched to + a T-vertex via a tight edge. """ # Assign label S to the blossom that contains vertex "v". - vb = matching.vertex_blossom[v] - assert stage_data.blossom_label[vb] == _LABEL_NONE - stage_data.blossom_label[vb] = _LABEL_S + bi = matching.vertex_blossom[i] + assert stage_data.blossom_label[bi] == _LABEL_NONE + stage_data.blossom_label[bi] = _LABEL_S - w = matching.vertex_mate[v] - if w == -1: - # Vertex "v" is unmatched. + j = matching.vertex_mate[i] + if j == -1: + # Vertex "i" is unmatched. # It must be either a top-level vertex or the base vertex of # a top-level blossom. - assert (vb == v) or (matching.get_blossom(vb).base_vertex == v) + assert (bi == i) or (matching.get_blossom(bi).base_vertex == i) # Mark the blossom that contains "v" as root of an alternating tree. - stage_data.blossom_link[vb] = None + stage_data.blossom_link[bi] = None else: - # Vertex "v" is matched to T-vertex "w". - wb = matching.vertex_blossom[w] - assert stage_data.blossom_label[wb] == _LABEL_T + # Vertex "i" is matched to T-vertex "j". + bj = matching.vertex_blossom[j] + assert stage_data.blossom_label[bj] == _LABEL_T - # Attach the blossom that contains "v" to the alternating tree. - stage_data.blossom_link[vb] = (w, v) + # Attach the blossom that contains "i" to the alternating tree. + stage_data.blossom_link[bi] = (j, i) # Initialize the list of least-slack edges of the newly labeled blossom. # This list will be filled by scanning the vertices of the blossom. - stage_data.blossom_best_edge_set[vb] = [] + stage_data.blossom_best_edge_set[bi] = [] # Add all vertices inside the newly labeled S-blossom to the queue. - if vb == v: - stage_data.queue.append(v) + if bi == i: + stage_data.queue.append(i) else: - stage_data.queue.extend(matching.blossom_vertices(vb)) + stage_data.queue.extend(matching.blossom_vertices(bi)) def _substage_assign_label_t( matching: _PartialMatching, stage_data: _StageData, - v: int, - w: int + i: int, + j: int ) -> None: - """Assign label T to the unlabeled blossom that contains vertex "w". + """Assign label T to the unlabeled blossom that contains vertex "j". - Attach it to the alternating tree via edge (v, w). - Then immediately assign label S to the mate of vertex "w". + Attach it to the alternating tree via edge (i, j). + Then immediately assign label S to the mate of vertex "j". Preconditions: - - "v" is an S-vertex. - - "w" is an unlabeled, matched vertex. - - There is a tight edge between vertices "v" and "w". + - "i" is an S-vertex. + - "j" is an unlabeled, matched vertex. + - There is a tight edge between vertices "i" and "j". """ - assert stage_data.blossom_label[matching.vertex_blossom[v]] == _LABEL_S + assert stage_data.blossom_label[matching.vertex_blossom[i]] == _LABEL_S # Assign label T to the unlabeled blossom. - wb = matching.vertex_blossom[w] - assert stage_data.blossom_label[wb] == _LABEL_NONE - stage_data.blossom_label[wb] = _LABEL_T - stage_data.blossom_link[wb] = (v, w) + bj = matching.vertex_blossom[j] + assert stage_data.blossom_label[bj] == _LABEL_NONE + stage_data.blossom_label[bj] = _LABEL_T + stage_data.blossom_link[bj] = (i, j) - # Assign label S to the blossom that contains the mate of vertex "w". - wbase = w if wb == w else matching.get_blossom(wb).base_vertex - y = matching.vertex_mate[wbase] - _substage_assign_label_s(matching, stage_data, y) + # Assign label S to the blossom that contains the mate of vertex "j". + jbase = j if bj == j else matching.get_blossom(bj).base_vertex + k = matching.vertex_mate[jbase] + assert k != -1 + _substage_assign_label_s(matching, stage_data, k) def _substage_add_s_to_s_edge( matching: _PartialMatching, stage_data: _StageData, - v: int, - w: int + i: int, + j: int ) -> Optional[_AlternatingPath]: - """Add the edge between S-vertices "v" and "w". + """Add the edge between S-vertices "i" and "j". If the edge connects blossoms that are part of the same alternating tree, this function creates a new S-blossom and returns None. @@ -933,10 +938,10 @@ def _substage_add_s_to_s_edge( Augmenting path if found; otherwise None. """ - # Trace back through the alternating trees from "v" and "w". - path = _trace_alternating_paths(matching, stage_data, v, w) + # Trace back through the alternating trees from "i" and "j". + path = _trace_alternating_paths(matching, stage_data, i, j) - # If the path is a cycle, we create a new blossom. + # If the path is a cycle, create a new blossom. # Otherwise the path is an augmenting path. # Note that an alternating starts and ends in the same blossom, # but not necessarily in the same vertex within that blossom. @@ -972,71 +977,70 @@ def _substage_scan( while stage_data.queue: # Take a vertex from the queue. - v = stage_data.queue.pop() + i = stage_data.queue.pop() # Double-check that "v" is an S-vertex. - vblossom = matching.vertex_blossom[v] - assert stage_data.blossom_label[vblossom] == _LABEL_S + bi = matching.vertex_blossom[i] + assert stage_data.blossom_label[bi] == _LABEL_S # Scan the edges that are incident on "v". - for e in adjacent_edges[v]: - (i, j, _wt) = edges[e] - w = i if i != v else j # TODO : consider abstracting this + for e in adjacent_edges[i]: + (p, q, _w) = edges[e] + j = p if p != i else q # TODO : consider abstracting this - # Consider the edge between vertices "v" and "w". + # Consider the edge between vertices "i" and "j". # Try to pull this edge into an alternating tree. - # Ignore edges that are internal to a blossom. - # Note: blossom index of vertex "v" may change during this loop, + # Note: blossom index of vertex "i" may change during this loop, # so we need to refresh it here. - vblossom = matching.vertex_blossom[v] - wblossom = matching.vertex_blossom[w] - if vblossom == wblossom: + bi = matching.vertex_blossom[i] + bj = matching.vertex_blossom[j] + + # Ignore edges that are internal to a blossom. + if bi == bj: continue - wlabel = stage_data.blossom_label[wblossom] + jlabel = stage_data.blossom_label[bj] # Check whether this edge is tight (has zero slack). # Only tight edges may be part of an alternating tree. slack = matching.edge_slack_2x(e) if slack <= 0: - if wlabel == _LABEL_NONE: - # Assign label T to the blossom that contains vertex "w". - _substage_assign_label_t(matching, stage_data, v, w) - elif wlabel == _LABEL_S: + if jlabel == _LABEL_NONE: + # Assign label T to the blossom that contains vertex "j". + _substage_assign_label_t(matching, stage_data, i, j) + elif jlabel == _LABEL_S: # This edge connects two S-blossoms. Use it to find # either a new blossom or an augmenting path. alternating_path = _substage_add_s_to_s_edge(matching, stage_data, - v, w) + i, j) if alternating_path is not None: return alternating_path - elif wlabel == _LABEL_S: + elif jlabel == _LABEL_S: # Update tracking of least-slack edges between S-blossoms. - best_edge = stage_data.blossom_best_edge[vblossom] + best_edge = stage_data.blossom_best_edge[bi] if ((best_edge < 0) or (slack < matching.edge_slack_2x(best_edge))): - stage_data.blossom_best_edge[vblossom] = e + stage_data.blossom_best_edge[bi] = e # Update the list of least-slack edges to S-blossoms for - # the blossom that contains "v". + # the blossom that contains "i". # We only do this for non-trivial blossoms. - if vblossom != v: - best_edge_set = stage_data.blossom_best_edge_set[vblossom] + if bi != i: + best_edge_set = stage_data.blossom_best_edge_set[bi] assert best_edge_set is not None best_edge_set.append(e) - if wlabel != _LABEL_S: - # Update tracking of least-slack edges from vertex "w" to + if jlabel != _LABEL_S: + # Update tracking of least-slack edges from vertex "j" to # any S-vertex. We do this only for T-vertices and unlabeled - # vertices. Edges which have zero slack are also tracked. - # If "w" is part of a T-blossom, it may become unlabeled - # later. At that point we will need a zero-slack edge to - # relabel vertex "w". - best_edge = stage_data.vertex_best_edge[w] + # vertices. Edges which already have zero slack are still + # tracked. + best_edge = stage_data.vertex_best_edge[j] if best_edge < 0 or slack < matching.edge_slack_2x(best_edge): - stage_data.vertex_best_edge[w] = e + stage_data.vertex_best_edge[j] = e # No further S vertices to scan, and no augmenting path found. return None @@ -1045,9 +1049,9 @@ def _substage_scan( def _find_path_through_blossom( matching: _PartialMatching, b: int, - s: int + sub: int ) -> tuple[list[int], list[tuple[int, int]]]: - """Construct a path through blossom "b" from sub-blossom "s" + """Construct a path through blossom "b" from sub-blossom "sub" to the base of the blossom. Return: @@ -1056,40 +1060,40 @@ def _find_path_through_blossom( blossom = matching.blossom[b] assert blossom is not None - nodes: list[int] = [s] + nodes: list[int] = [sub] edges: list[tuple[int, int]] = [] # Walk around the blossom from "s" to its base. - i = blossom.subblossoms.index(s) + p = blossom.subblossoms.index(sub) nsub = len(blossom.subblossoms) - while i != 0: - if i % 2 == 0: + while p != 0: + if p % 2 == 0: # Stepping towards the beginning of the subblossom list. - # Currently at subblossom (i), next position (i-2): + # Currently at subblossom (p), next position (p-2): # - # 0 --- 1 === 2 --- 3 === (i-2) --- (i-1) ==(v,w)== (i) + # 0 --- 1 === 2 --- 3 === (p-2) --- (p-1) ==(i,j)== (p) # ^^^ ^^^ # <------------------- # - # We flip edges from (v,w) to (w,v) to make them fit + # We flip edges from (i,j) to (j,i) to make them fit # in the path from "s" to base. - edges.append(blossom.edges[i-1][::-1]) - nodes.append(blossom.subblossoms[i-1]) - edges.append(blossom.edges[i-2][::-1]) - nodes.append(blossom.subblossoms[i-2]) - i -= 2 + edges.append(blossom.edges[p-1][::-1]) + nodes.append(blossom.subblossoms[p-1]) + edges.append(blossom.edges[p-2][::-1]) + nodes.append(blossom.subblossoms[p-2]) + p -= 2 else: # Stepping towards the end of the subblossom list. - # Currently at subblossom (i), next position (i+2): + # Currently at subblossom (p), next position (p+2): # - # (i) ==(v,w)== (i+1) --- (i+2) === (i+3) --- 0 + # (p) ==(i,j)== (p+1) --- (p+2) === (p+3) --- 0 # ^^^ ^^^ # -------------------> - edges.append(blossom.edges[i]) - nodes.append(blossom.subblossoms[i+1]) - edges.append(blossom.edges[i+1]) - nodes.append(blossom.subblossoms[(i+2) % nsub]) - i = (i + 2) % nsub + edges.append(blossom.edges[p]) + nodes.append(blossom.subblossoms[p+1]) + edges.append(blossom.edges[p+1]) + nodes.append(blossom.subblossoms[(p+2) % nsub]) + p = (p + 2) % nsub return (nodes, edges) @@ -1117,8 +1121,8 @@ def _expand_t_blossom( if sub < num_vertex: matching.vertex_blossom[sub] = sub else: - for v in matching.blossom_vertices(sub): - matching.vertex_blossom[v] = sub + for i in matching.blossom_vertices(sub): + matching.vertex_blossom[i] = sub assert stage_data.blossom_label[sub] == _LABEL_NONE # The expanding blossom was part of an alternating tree, linked to @@ -1127,38 +1131,39 @@ def _expand_t_blossom( # We must reconstruct this part of the alternating tree, which will # now run via sub-blossoms of the expanded blossom. - # Determine which sub-blossom is attached to the parent tree node. + # Find the sub-blossom that is attached to the parent node in + # the alternating tree. # TODO : uglyness with the assertion entry_link = stage_data.blossom_link[b] assert entry_link is not None - (v, w) = entry_link - sub = matching.vertex_blossom[w] + (i, j) = entry_link + sub = matching.vertex_blossom[j] # Assign label T to that sub-blossom. stage_data.blossom_label[sub] = _LABEL_T - stage_data.blossom_link[sub] = (v, w) + stage_data.blossom_link[sub] = (i, j) # Walk through the expanded blossom from "sub" to the base vertex. # Assign alternating S and T labels to the sub-blossoms and attach # them to the alternating tree. (path_nodes, path_edges) = _find_path_through_blossom(matching, b, sub) - for i in range(0, len(path_edges), 2): + for p in range(0, len(path_edges), 2): # - # (i) ===== (i+1) ----- (i+2) - # T S T + # (p) ==(j,i)== (p+1) ----- (p+2) + # T S T # - # path_nodes[i] has already been labeled T. - # We now assign labels to path_nodes[i+1] and path_nodes[i+2]. + # path_nodes[p] has already been labeled T. + # We now assign labels to path_nodes[p+1] and path_nodes[p+2]. - # Assign label S to path_nodes[i+1]. - (w, v) = path_edges[i] - _substage_assign_label_s(matching, stage_data, v) + # Assign label S to path_nodes[p+1]. + (j, i) = path_edges[p] + _substage_assign_label_s(matching, stage_data, i) - # Assign label T to path_nodes[i+2] and attach it to path_nodes[i+1]. - sub = path_nodes[i+2] + # Assign label T to path_nodes[i+2] and attach it to path_nodes[p+1]. + sub = path_nodes[p+2] stage_data.blossom_label[sub] = _LABEL_T - stage_data.blossom_link[sub] = path_edges[i+1] + stage_data.blossom_link[sub] = path_edges[p+1] # Unlabel and delete the expanded blossom. Recycle its blossom index. stage_data.blossom_label[b] = _LABEL_NONE @@ -1208,22 +1213,23 @@ def _expand_zero_dual_blossoms(matching: _PartialMatching) -> None: matching.vertex_blossom[sub] = sub else: # Non-trivial sub-blossom. - # If its dual variable is zero, we expand it recursively. + # If its dual variable is zero, expand it recursively. if matching.get_blossom(sub).dual_var == 0: stack.append(sub) else: - # This sub-blossom will not be expanded. - # We still need to update its "vertex_blossom" entries. - for v in matching.blossom_vertices(sub): - matching.vertex_blossom[v] = sub + # This sub-blossom will not be expanded; + # it now becomes top-level. Update its vertices + # to point to this sub-blossom. + for i in matching.blossom_vertices(sub): + matching.vertex_blossom[i] = sub # Delete the expanded blossom. Recycle its blossom index. matching.blossom[b] = None matching.unused_blossoms.append(b) -def _augment_blossom(matching: _PartialMatching, b: int, v: int) -> None: - """Augment along an alternating path through blossom "b", from vertex "v" +def _augment_blossom(matching: _PartialMatching, b: int, i: int) -> None: + """Augment along an alternating path through blossom "b", from vertex "i" to the base vertex of the blossom. This function takes time O(n). @@ -1231,7 +1237,7 @@ def _augment_blossom(matching: _PartialMatching, b: int, v: int) -> None: num_vertex = matching.graph.num_vertex # Use an explicit stack to avoid deep recursion. - stack = [(b, v)] + stack = [(b, i)] while stack: (top_blossom, sub) = stack.pop() @@ -1248,41 +1254,41 @@ def _augment_blossom(matching: _PartialMatching, b: int, v: int) -> None: assert blossom is not None # Walk through the expanded blossom from "sub" to the base vertex. - (path_nodes, path_edges) = _find_path_through_blossom( - matching, b, sub) + (path_nodes, path_edges) = _find_path_through_blossom(matching, + b, sub) - for i in range(0, len(path_edges), 2): + for p in range(0, len(path_edges), 2): # Before augmentation: - # path_nodes[i] is matched to path_nodes[i+1] + # path_nodes[p] is matched to path_nodes[p+1] # - # (i) ===== (i+1) ---(v,w)--- (i+2) + # (p) ===== (p+1) ---(i,j)--- (p+2) # # After augmentation: - # path_nodes[i+1] is mathed to path_nodes[i+2] via edge (v, w) + # path_nodes[p+1] is mathed to path_nodes[p+2] via edge (i, j) # - # (i) ----- (i+1) ===(v,w)=== (i+2) + # (p) ----- (p+1) ===(i,j)=== (p+2) # - # Pull the edge (v, w) into the matching. - (v, w) = path_edges[i+1] - matching.vertex_mate[v] = w - matching.vertex_mate[w] = v + # Pull the edge (i, j) into the matching. + (i, j) = path_edges[p+1] + matching.vertex_mate[i] = j + matching.vertex_mate[j] = i - # Augment through the subblossoms touching the edge (v, w). + # Augment through the subblossoms touching the edge (i, j). # Nothing needs to be done for trivial subblossoms. - vb = path_nodes[i+1] - if vb >= num_vertex: - stack.append((vb, v)) + bi = path_nodes[p+1] + if bi >= num_vertex: + stack.append((bi, i)) - wb = path_nodes[i+2] - if wb >= num_vertex: - stack.append((wb, w)) + bj = path_nodes[p+2] + if bj >= num_vertex: + stack.append((bj, j)) # Rotate the subblossom list so the new base ends up in position 0. - sub_pos = blossom.subblossoms.index(sub) - blossom.subblossoms = (blossom.subblossoms[sub_pos:] - + blossom.subblossoms[:sub_pos]) - blossom.edges = blossom.edges[sub_pos:] + blossom.edges[:sub_pos] + p = blossom.subblossoms.index(sub) + blossom.subblossoms = ( + blossom.subblossoms[p:] + blossom.subblossoms[:p]) + blossom.edges = blossom.edges[p:] + blossom.edges[:p] # Update the base vertex. # We can pull this from the sub-blossom where we started since @@ -1305,28 +1311,28 @@ def _augment_matching( # Check that the augmenting path starts and ends in an unmatched vertex # or a blossom with unmatched base. assert len(path.edges) % 2 == 1 - for v in (path.edges[0][0], path.edges[-1][1]): - b = matching.vertex_blossom[v] - if b != v: - v = matching.get_blossom(b).base_vertex - assert matching.vertex_mate[v] == -1 + for i in (path.edges[0][0], path.edges[-1][1]): + b = matching.vertex_blossom[i] + if b != i: + i = matching.get_blossom(b).base_vertex + assert matching.vertex_mate[i] == -1 # Consider the edges of the augmenting path that are currently not # part of the matching but will become part of it. - for (v, w) in path.edges[0::2]: + for (i, j) in path.edges[0::2]: # Augment through non-trivial blossoms on either side of this edge. - vblossom = matching.vertex_blossom[v] - if vblossom != v: - _augment_blossom(matching, vblossom, v) + bi = matching.vertex_blossom[i] + if bi != i: + _augment_blossom(matching, bi, i) - wblossom = matching.vertex_blossom[w] - if wblossom != w: - _augment_blossom(matching, wblossom, w) + bj = matching.vertex_blossom[j] + if bj != j: + _augment_blossom(matching, bj, j) # Pull the edge into the matching. - matching.vertex_mate[v] = w - matching.vertex_mate[w] = v + matching.vertex_mate[i] = j + matching.vertex_mate[j] = i def _calc_dual_delta( @@ -1356,16 +1362,17 @@ def _calc_dual_delta( # Compute delta1: minimum dual variable of any S-vertex. delta_type = 1 - delta_2x = min(matching.dual_var_2x[v] - for v in range(num_vertex) - if stage_data.blossom_label[matching.vertex_blossom[v]]) + delta_2x = min( + matching.dual_var_2x[i] + for i in range(num_vertex) + if stage_data.blossom_label[matching.vertex_blossom[i]] == _LABEL_S) # Compute delta2: minimum slack of any edge between an S-vertex and # an unlabeled vertex. - for v in range(num_vertex): - vb = matching.vertex_blossom[v] - if stage_data.blossom_label[vb] == _LABEL_NONE: - e = stage_data.vertex_best_edge[v] + for i in range(num_vertex): + bi = matching.vertex_blossom[i] + if stage_data.blossom_label[bi] == _LABEL_NONE: + e = stage_data.vertex_best_edge[i] if e != -1: slack = matching.edge_slack_2x(e) if slack <= delta_2x: @@ -1417,14 +1424,14 @@ def _apply_delta_step( num_vertex = matching.graph.num_vertex # Apply delta to dual variables of all vertices. - for v in range(num_vertex): - vlabel = stage_data.blossom_label[matching.vertex_blossom[v]] - if vlabel == _LABEL_S: + for i in range(num_vertex): + ilabel = stage_data.blossom_label[matching.vertex_blossom[i]] + if ilabel == _LABEL_S: # S-vertex: subtract delta from dual variable. - matching.dual_var_2x[v] -= delta_2x - elif vlabel == _LABEL_T: + matching.dual_var_2x[i] -= delta_2x + elif ilabel == _LABEL_T: # T-vertex: add delta to dual variable. - matching.dual_var_2x[v] += delta_2x + matching.dual_var_2x[i] += delta_2x # Apply delta to dual variables of top-level non-trivial blossoms. for b in range(num_vertex, 2 * num_vertex): @@ -1459,9 +1466,9 @@ def _run_stage(matching: _PartialMatching) -> bool: stage_data = _StageData(matching.graph) # Assign label S to all unmatched vertices and put them in the queue. - for v in range(num_vertex): - if matching.vertex_mate[v] == -1: - _substage_assign_label_s(matching, stage_data, v) + for i in range(num_vertex): + if matching.vertex_mate[i] == -1: + _substage_assign_label_s(matching, stage_data, i) # Stop if all vertices are matched. # No further improvement is possible in this case. @@ -1494,18 +1501,18 @@ def _run_stage(matching: _PartialMatching) -> bool: if delta_type == 2: # Use the edge from S-vertex to unlabeled vertex that got # unlocked through the delta update. - (v, w, _wt) = matching.graph.edges[delta_edge] - if (stage_data.blossom_label[matching.vertex_blossom[v]] != + (i, j, _w) = matching.graph.edges[delta_edge] + if (stage_data.blossom_label[matching.vertex_blossom[i]] != _LABEL_S): - (v, w) = (w, v) - _substage_assign_label_t(matching, stage_data, v, w) + (i, j) = (j, i) + _substage_assign_label_t(matching, stage_data, i, j) elif delta_type == 3: # Use the S-to-S edge that got unlocked through the delta update. # This may reveal an augmenting path. - (v, w, _wt) = matching.graph.edges[delta_edge] + (i, j, _w) = matching.graph.edges[delta_edge] augmenting_path = _substage_add_s_to_s_edge( - matching, stage_data, v, w) + matching, stage_data, i, j) if augmenting_path is not None: break @@ -1553,12 +1560,12 @@ def _verify_optimum(matching: _PartialMatching) -> None: # Double-check that each matching edge actually exists in the graph. num_matched_vertex = 0 - for v in range(num_vertex): - if vertex_mate[v] != -1: + for i in range(num_vertex): + if vertex_mate[i] != -1: num_matched_vertex += 1 num_matched_edge = 0 - for (i, j, _wt) in matching.graph.edges: + for (i, j, _w) in matching.graph.edges: if vertex_mate[i] == j: num_matched_edge += 1 @@ -1570,8 +1577,8 @@ def _verify_optimum(matching: _PartialMatching) -> None: # Count the number of vertices in each blossom. blossom_nvertex = (2 * num_vertex) * [0] - for v in range(num_vertex): - b = matching.blossom_parent[v] + for i in range(num_vertex): + b = matching.blossom_parent[i] while b != -1: blossom_nvertex[b] += 1 b = matching.blossom_parent[b] @@ -1580,7 +1587,7 @@ def _verify_optimum(matching: _PartialMatching) -> None: # Also count the number of matched edges in each blossom. blossom_nmatched = (2 * num_vertex) * [0] - for (i, j, wt) in matching.graph.edges: + for (i, j, w) in matching.graph.edges: # List blossoms that contain vertex "i". iblossoms = [] @@ -1608,7 +1615,7 @@ def _verify_optimum(matching: _PartialMatching) -> None: # + sum(dual[b] for blossoms "b" containing the edge) # # Multiply weights by 2 to ensure integer values. - slack = vertex_dual_var_2x[i] + vertex_dual_var_2x[j] - 2 * wt + slack = vertex_dual_var_2x[i] + vertex_dual_var_2x[j] - 2 * w slack += 2 * sum(blossom_dual_var[b] for b in edge_blossoms) # Check that all edges have non-negative slack. @@ -1624,9 +1631,9 @@ def _verify_optimum(matching: _PartialMatching) -> None: blossom_nmatched[b] += 1 # Check that all unmatched vertices have zero dual. - for v in range(num_vertex): - if vertex_mate[v] == -1: - assert vertex_dual_var_2x[v] == 0 + for i in range(num_vertex): + if vertex_mate[i] == -1: + assert vertex_dual_var_2x[i] == 0 # Check that all blossoms with positive dual are "full". # A blossom is full if all except one of its vertices are matched