1
0
Fork 0
maximum-weight-matching/python/mwmatching.py

2117 lines
79 KiB
Python
Raw Normal View History

2023-02-04 21:06:34 +01:00
"""
Algorithm for finding a maximum weight matching in general graphs.
"""
from __future__ import annotations
import sys
2024-06-25 21:01:52 +02:00
import itertools
2023-02-04 21:06:34 +01:00
import math
from collections.abc import Sequence
from typing import NamedTuple, Optional
2023-02-04 21:06:34 +01:00
from datastruct import UnionFindQueue, PriorityQueue
2023-02-04 21:06:34 +01:00
def maximum_weight_matching(
edges: Sequence[tuple[int, int, float]]
2023-02-04 21:06:34 +01:00
) -> list[tuple[int, int]]:
"""Compute a maximum-weighted matching in the general undirected weighted
graph given by "edges".
The graph is specified as a list of edges, each edge specified as a tuple
of its two vertices and the edge weight.
There may be at most one edge between any pair of vertices.
No vertex may have an edge to itself.
The graph may be non-connected (i.e. contain multiple components).
Vertices are indexed by consecutive, non-negative integers, such that
the first vertex has index 0 and the last vertex has index (n-1).
Edge weights may be integers or floating point numbers.
Isolated vertices (not incident to any edge) are allowed, but not
recommended since such vertices consume time and memory but have
no effect on the maximum-weight matching.
Edges with negative weight are ignored.
2023-02-04 21:06:34 +01:00
This function takes time O(n * (n + m) * log(n)),
where "n" is the number of vertices and "m" is the number of edges.
This function uses O(n + m) memory.
2023-02-04 21:06:34 +01:00
Parameters:
2023-02-07 11:23:42 +01:00
edges: List of edges, each edge specified as a tuple "(x, y, w)"
where "x" and "y" are vertex indices and "w" is the edge weight.
2023-02-04 21:06:34 +01:00
Returns:
List of pairs of matched vertex indices.
2023-02-07 11:23:42 +01:00
This is a subset of the edges in the graph.
It contains a tuple "(x, y)" if vertex "x" is matched to vertex "y".
2023-02-04 21:06:34 +01:00
Raises:
ValueError: If the input does not satisfy the constraints.
TypeError: If the input contains invalid data types.
2023-02-12 21:18:57 +01:00
MatchingError: If the matching algorithm fails.
This can only happen if there is a bug in the algorithm.
2023-02-04 21:06:34 +01:00
"""
# Check that the input meets all constraints.
_check_input_types(edges)
_check_input_graph(edges)
2023-02-06 19:58:27 +01:00
# Remove edges with negative weight.
edges = _remove_negative_weight_edges(edges)
2023-02-04 21:06:34 +01:00
# Special case for empty graphs.
if not edges:
return []
# Initialize graph representation.
graph = _GraphInfo(edges)
2023-02-06 22:33:26 +01:00
# Initialize the matching algorithm.
ctx = _MatchingContext(graph)
ctx.start()
2023-02-04 21:06:34 +01:00
# Improve the solution until no further improvement is possible.
#
# Each successful pass through this loop increases the number
# of matched edges by 1.
#
# This loop runs through at most (n/2 + 1) iterations.
2024-06-29 22:06:28 +02:00
# Each iteration takes time O((n + m) * log(n)).
2023-02-06 22:33:26 +01:00
while ctx.run_stage():
2023-02-04 21:06:34 +01:00
pass
# Extract the final solution.
ctx.cleanup()
2023-02-04 21:06:34 +01:00
pairs: list[tuple[int, int]] = [
2023-02-07 11:23:42 +01:00
(x, y) for (x, y, _w) in edges if ctx.vertex_mate[x] == y]
2023-02-04 21:06:34 +01:00
# Verify that the matching is optimal.
# This is just a safeguard; the verification will always pass unless
# there is a bug in the matching algorithm.
# Verification only works reliably for integer weights.
2023-02-04 21:06:34 +01:00
if graph.integer_weights:
2023-02-06 22:33:26 +01:00
_verify_optimum(ctx)
2023-02-04 21:06:34 +01:00
return pairs
def adjust_weights_for_maximum_cardinality_matching(
edges: Sequence[tuple[int, int, float]]
) -> Sequence[tuple[int, int, float]]:
2023-02-04 21:06:34 +01:00
"""Adjust edge weights such that the maximum-weight matching of
the adjusted graph is a maximum-cardinality matching, equal to
a matching in the original graph that has maximum weight out of all
matchings with maximum cardinality.
The graph is specified as a list of edges, each edge specified as a tuple
of its two vertices and the edge weight.
Edge weights may be integers or floating point numbers.
Negative edge weights are allowed.
This function increases all edge weights by an equal amount such that
the adjusted weights satisfy the following conditions:
- All edge weights are positive;
- The minimum edge weight is at least "n" times the difference between
maximum and minimum edge weight.
2023-02-13 22:38:40 +01:00
These conditions ensure that a maximum-cardinality matching will be found.
Proof: The weight of any non-maximum-cardinality matching can be increased
by matching an additional edge, even if the new edge has minimum edge
weight and causes all other matched edges to degrade from maximum to
minimum edge weight.
2023-02-04 21:06:34 +01:00
Since we are only considering maximum-cardinality matchings, increasing
all edge weights by an equal amount will not change the set of edges
that makes up the maximum-weight matching.
This function increases edge weights by an amount that is proportional
to the product of the unadjusted weight range and the number of vertices
in the graph. In case of a big graph with floating point weights, this
may introduce rounding errors in the weights.
This function takes time O(m), where "m" is the number of edges.
Parameters:
2023-02-07 11:23:42 +01:00
edges: List of edges, each edge specified as a tuple "(x, y, w)"
where "x" and "y" are vertex indices and "w" is the edge weight.
2023-02-04 21:06:34 +01:00
Returns:
List of edges with adjusted weights. If no adjustments are necessary,
the input list instance may be returned.
Raises:
ValueError: If the input does not satisfy the constraints.
TypeError: If the input contains invalid data types.
"""
_check_input_types(edges)
# Don't worry about empty graphs:
if not edges:
return edges
2023-02-07 11:23:42 +01:00
num_vertex = 1 + max(max(x, y) for (x, y, _w) in edges)
2023-02-04 21:06:34 +01:00
2023-02-07 11:23:42 +01:00
min_weight = min(w for (_x, _y, w) in edges)
max_weight = max(w for (_x, _y, w) in edges)
2023-02-04 21:06:34 +01:00
weight_range = max_weight - min_weight
# Do nothing if the weights already ensure a maximum-cardinality matching.
if min_weight > 0 and min_weight >= num_vertex * weight_range:
return edges
delta: float
2023-02-04 21:06:34 +01:00
if weight_range > 0:
# Increase weights to make minimum edge weight large enough
# to improve any non-maximum-cardinality matching.
delta = num_vertex * weight_range - min_weight
else:
# All weights are the same. Increase weights to make them positive.
delta = 1 - min_weight
assert delta >= 0
# Increase all edge weights by "delta".
2023-02-07 11:23:42 +01:00
return [(x, y, w + delta) for (x, y, w) in edges]
2023-02-04 21:06:34 +01:00
2023-02-12 21:18:57 +01:00
class MatchingError(Exception):
"""Raised when verification of the matching fails.
This can only happen if there is a bug in the algorithm.
"""
pass
def _check_input_types(edges: Sequence[tuple[int, int, float]]) -> None:
2023-02-04 21:06:34 +01:00
"""Check that the input consists of valid data types and valid
numerical ranges.
This function takes time O(m).
Parameters:
2023-02-07 11:23:42 +01:00
edges: List of edges, each edge specified as a tuple "(x, y, w)"
where "x" and "y" are edge indices and "w" is the edge weight.
2023-02-04 21:06:34 +01:00
Raises:
ValueError: If the input does not satisfy the constraints.
TypeError: If the input contains invalid data types.
"""
float_limit = sys.float_info.max / 4
if not isinstance(edges, list):
raise TypeError('"edges" must be a list')
for e in edges:
if (not isinstance(e, tuple)) or (len(e) != 3):
raise TypeError("Each edge must be specified as a 3-tuple")
2023-02-07 11:23:42 +01:00
(x, y, w) = e
2023-02-04 21:06:34 +01:00
2023-02-07 11:23:42 +01:00
if (not isinstance(x, int)) or (not isinstance(y, int)):
2023-02-04 21:06:34 +01:00
raise TypeError("Edge endpoints must be integers")
2023-02-07 11:23:42 +01:00
if (x < 0) or (y < 0):
2023-02-04 21:06:34 +01:00
raise ValueError("Edge endpoints must be non-negative integers")
2023-02-06 17:39:03 +01:00
if not isinstance(w, (int, float)):
2023-02-04 21:06:34 +01:00
raise TypeError(
"Edge weights must be integers or floating point numbers")
2023-02-06 17:39:03 +01:00
if isinstance(w, float):
if not math.isfinite(w):
2023-02-04 21:06:34 +01:00
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.
2023-02-06 17:39:03 +01:00
if w > float_limit:
2023-02-04 21:06:34 +01:00
raise ValueError("Floating point edge weights must be"
f" less than {float_limit:g}")
def _check_input_graph(edges: Sequence[tuple[int, int, float]]) -> None:
2023-02-04 21:06:34 +01:00
"""Check that the input is a valid graph, without any multi-edges and
without any self-edges.
This function takes time O(m * log(m)).
Parameters:
2023-02-07 11:23:42 +01:00
edges: List of edges, each edge specified as a tuple "(x, y, w)"
where "x" and "y" are edge indices and "w" is the edge weight.
2023-02-04 21:06:34 +01:00
Raises:
ValueError: If the input does not satisfy the constraints.
"""
# Check that the graph has no self-edges.
2023-02-07 11:23:42 +01:00
for (x, y, _w) in edges:
if x == y:
2023-02-04 21:06:34 +01:00
raise ValueError("Self-edges are not supported")
# Check that the graph does not have multi-edges.
# 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.
2023-02-07 11:23:42 +01:00
edge_endpoints = [((x, y) if (x < y) else (y, x)) for (x, y, _w) in edges]
2023-02-04 21:06:34 +01:00
edge_endpoints.sort()
for i in range(len(edge_endpoints) - 1):
if edge_endpoints[i] == edge_endpoints[i+1]:
raise ValueError(f"Duplicate edge {edge_endpoints[i]}")
2023-02-06 19:58:27 +01:00
def _remove_negative_weight_edges(
edges: Sequence[tuple[int, int, float]]
) -> Sequence[tuple[int, int, float]]:
2023-02-06 19:58:27 +01:00
"""Remove edges with negative weight.
This does not change the solution of the maximum-weight matching problem,
but prevents complications in the algorithm.
"""
if any(e[2] < 0 for e in edges):
return [e for e in edges if e[2] >= 0]
else:
return edges
2023-02-04 21:06:34 +01:00
class _GraphInfo:
"""Representation of the input graph.
These data remain unchanged while the algorithm runs.
"""
def __init__(self, edges: Sequence[tuple[int, int, float]]) -> None:
2023-02-04 21:06:34 +01:00
"""Initialize the graph representation and prepare an adjacency list.
This function takes time O(n + m).
"""
# Vertices are indexed by integers in range 0 .. n-1.
# Edges are indexed by integers in range 0 .. m-1.
#
# Each edge is incident on two vertices.
# Each edge also has a weight.
#
2023-02-07 11:23:42 +01:00
# "edges[e] = (x, y, w)" where
2023-02-04 21:06:34 +01:00
# "e" is an edge index;
2023-02-07 11:23:42 +01:00
# "x" and "y" are vertex indices of the incident vertices;
2023-02-06 17:39:03 +01:00
# "w" is the edge weight.
2023-02-04 21:06:34 +01:00
#
# These data remain unchanged while the algorithm runs.
self.edges: Sequence[tuple[int, int, float]] = edges
2023-02-04 21:06:34 +01:00
# num_vertex = the number of vertices.
if edges:
2023-02-07 11:23:42 +01:00
self.num_vertex = 1 + max(max(x, y) for (x, y, _w) in edges)
2023-02-04 21:06:34 +01:00
else:
self.num_vertex = 0
# Each vertex is incident to zero or more edges.
#
# "adjacent_edges[x]" is the list of edge indices of edges incident
# to the vertex with index "x".
2023-02-04 21:06:34 +01:00
#
# These data remain unchanged while the algorithm runs.
self.adjacent_edges: list[list[int]] = [
[] for v in range(self.num_vertex)]
2023-02-07 11:23:42 +01:00
for (e, (x, y, _w)) in enumerate(edges):
self.adjacent_edges[x].append(e)
self.adjacent_edges[y].append(e)
2023-02-04 21:06:34 +01:00
# Determine whether _all_ weights are integers.
# In this case we can avoid floating point computations entirely.
2023-02-06 17:39:03 +01:00
self.integer_weights: bool = all(isinstance(w, int)
2023-02-07 11:23:42 +01:00
for (_x, _y, w) in edges)
2023-02-04 21:06:34 +01:00
2023-02-06 22:33:26 +01:00
# Each vertex may be labeled "S" (outer) or "T" (inner) or be unlabeled.
_LABEL_NONE = 0
_LABEL_S = 1
_LABEL_T = 2
2023-02-04 21:06:34 +01:00
class _Blossom:
"""Represents a blossom in a partially matched graph.
2023-02-04 21:06:34 +01:00
A blossom is an odd-length alternating cycle over sub-blossoms.
An alternating path consists of alternating matched and unmatched edges.
An alternating cycle is an alternating path that starts and ends in
the same sub-blossom.
A single vertex by itself is also a blossom: a "trivial blossom".
An instance of this class represents either a trivial blossom,
or a non-trivial blossom which consists of multiple sub-blossoms.
2023-02-04 21:06:34 +01:00
Blossoms are recursive structures: A non-trivial blossoms contains
sub-blossoms, which may themselves contain sub-blossoms etc.
2023-02-05 19:33:13 +01:00
Each blossom contains exactly one vertex that is not matched to another
vertex in the same blossom. This is the "base vertex" of the blossom.
"""
def __init__(self, base_vertex: int) -> None:
"""Initialize a new blossom."""
# If this is not a top-level blossom,
# "parent" is the blossom in which this blossom is a sub-blossom.
#
# If this is a top-level blossom,
# "parent = None".
self.parent: Optional[_NonTrivialBlossom] = None
# "base_vertex" is the vertex index of the base of the blossom.
# This is the unique vertex which is contained in the blossom
# but not matched to another vertex in the same blossom.
#
# For trivial blossoms, the base vertex is simply the only
# vertex in the blossom.
self.base_vertex: int = base_vertex
2024-06-30 21:38:12 +02:00
# A top-level blossom that is part of an alternating tree,
# has label S or T. An unlabeled top-level blossom is not part
# of any alternating tree.
self.label: int = _LABEL_NONE
2024-06-30 21:38:12 +02:00
# A labeled top-level blossoms keeps track of the edge through which
# it is attached to the alternating tree.
#
# "tree_edge = (x, y)" if the blossom is attached to an alternating
# tree via edge "(x, y)" and vertex "y" is contained in the blossom.
#
# "tree_edge = None" if the blossom is the root of an alternating tree.
self.tree_edge: Optional[tuple[int, int]] = None
# For a labeled top-level blossom,
2024-06-30 09:11:58 +02:00
# "tree_blossoms" is the set of all top-level blossoms that belong
# to the same alternating tree. The same set instance is shared by
# all top-level blossoms in the tree.
self.tree_blossoms: "Optional[set[_Blossom]]" = None
# Each top-level blossom maintains a union-find datastructure
# containing all vertices in the blossom.
self.vertex_set: "UnionFindQueue[_Blossom, int]"
self.vertex_set = UnionFindQueue(self)
# If this is a top-level unlabeled blossom with an edge to an
# S-blossom, "delta2_node" is the corresponding node in the delta2
# queue.
self.delta2_node: Optional[PriorityQueue.Node] = None
2024-06-30 21:38:12 +02:00
# This variable holds pending lazy updates to the dual variables
# of the vertices inside the blossom.
self.vertex_dual_offset: float = 0
# "marker" is a temporary variable used to discover common
2024-06-25 21:01:52 +02:00
# ancestors in the alternating tree. It is normally False, except
# when used by "trace_alternating_paths()".
self.marker: bool = False
def vertices(self) -> list[int]:
"""Return a list of vertex indices contained in the blossom."""
return [self.base_vertex]
class _NonTrivialBlossom(_Blossom):
"""Represents a non-trivial blossom in a partially matched graph.
A non-trivial blossom is a blossom that contains multiple sub-blossoms
(at least 3 sub-blossoms, since all blossoms have odd length).
Non-trivial blossoms maintain a list of their sub-blossoms and the edges
between their subblossoms.
Unlike trivial blossoms, each non-trivial blossom is associated with
a variable in the dual LPP problem.
Non-trivial blossoms are created and destroyed by the matching algorithm.
2023-02-04 21:06:34 +01:00
This implies that not every odd-length alternating cycle is a blossom;
it only becomes a blossom through an explicit action of the algorithm.
An existing blossom may change when the matching is augmented along
a path that runs through the blossom.
2023-02-04 21:06:34 +01:00
"""
def __init__(
self,
subblossoms: list[_Blossom],
edges: list[tuple[int, int]]
2023-02-04 21:06:34 +01:00
) -> None:
"""Initialize a new blossom."""
super().__init__(subblossoms[0].base_vertex)
2023-02-04 21:06:34 +01:00
# Sanity check.
n = len(subblossoms)
assert len(edges) == n
assert n >= 3
assert n % 2 == 1
# "subblossoms" is a list of the sub-blossoms of the blossom,
# ordered by their appearance in the alternating cycle.
#
# "subblossoms[0]" is the start and end of the alternating cycle.
# "subblossoms[0]" contains the base vertex of the blossom.
self.subblossoms: list[_Blossom] = subblossoms
2023-02-04 21:06:34 +01:00
# "edges" is a list of edges linking the sub-blossoms.
2023-02-07 11:23:42 +01:00
# Each edge is represented as an ordered pair "(x, y)" where "x"
# and "y" are vertex indices.
2023-02-04 21:06:34 +01:00
#
2023-02-07 11:23:42 +01:00
# "edges[0] = (x, y)" where vertex "x" in "subblossoms[0]" is
# adjacent to vertex "y" in "subblossoms[1]", etc.
2023-02-04 21:06:34 +01:00
self.edges: list[tuple[int, int]] = edges
# Every non-trivial blossom has a variable in the dual LPP.
2024-05-26 20:55:54 +02:00
# New blossoms start with dual variable 0.
2023-02-04 21:06:34 +01:00
#
2024-05-26 20:55:54 +02:00
# The value of the dual variable changes through delta steps,
# but these changes are implemented as lazy updates.
2024-05-26 00:07:16 +02:00
#
2024-05-26 20:55:54 +02:00
# The true dual value of a top-level S-blossom is
# blossom.dual_var + ctx.delta_sum_2x
2024-05-26 11:32:02 +02:00
#
2024-05-26 20:55:54 +02:00
# The true dual value of a top-level T-blossom is
# blossom.dual_var - ctx.delta_sum_2x
2024-05-26 00:07:16 +02:00
#
2024-05-26 20:55:54 +02:00
# The true dual value of any other type of blossom is simply
# blossom.dual_var
2024-05-26 11:32:02 +02:00
#
2024-05-26 20:55:54 +02:00
# Note that "dual_var" is invariant under delta steps.
self.dual_var: float = 0
2023-02-04 21:06:34 +01:00
# If this is a top-level T-blossom,
# "delta4_node" is the corresponding node in the delta4 queue.
# Otherwise "delta4_node" is None.
self.delta4_node: Optional[PriorityQueue.Node] = None
def vertices(self) -> list[int]:
"""Return a list of vertex indices contained in the blossom."""
# Use an explicit stack to avoid deep recursion.
stack: list[_NonTrivialBlossom] = [self]
nodes: list[int] = []
while stack:
b = stack.pop()
for sub in b.subblossoms:
if isinstance(sub, _NonTrivialBlossom):
stack.append(sub)
else:
nodes.append(sub.base_vertex)
return nodes
2023-02-04 21:06:34 +01:00
2023-02-06 22:33:26 +01:00
class _AlternatingPath(NamedTuple):
"""Represents a list of edges forming an alternating path or an
alternating cycle."""
2023-02-06 22:33:26 +01:00
edges: list[tuple[int, int]]
is_cycle: bool
2023-02-06 22:33:26 +01:00
class _MatchingContext:
"""Holds all data used by the matching algorithm.
It contains a partial solution of the matching problem and several
auxiliary data structures.
"""
def __init__(self, graph: _GraphInfo) -> None:
"""Set up the initial state of the matching algorithm."""
num_vertex = graph.num_vertex
# Reference to the input graph.
# The graph does not change while the algorithm runs.
self.graph = graph
# Each vertex is either single (unmatched) or matched to
# another vertex.
#
# If vertex "x" is matched to vertex "y",
# "vertex_mate[x] == y" and "vertex_mate[y] == x".
#
# If vertex "x" is unmatched, "vertex_mate[x] == -1".
#
# Initially all vertices are unmatched.
self.vertex_mate: list[int] = num_vertex * [-1]
# Each vertex is associated with a trivial blossom.
# In addition, non-trivial blossoms may be created and destroyed
# during the course of the matching algorithm.
#
# "trivial_blossom[x]" is the trivial blossom that contains only
# vertex "x".
self.trivial_blossom: list[_Blossom] = [_Blossom(x)
for x in range(num_vertex)]
# Non-trivial blossoms may be created and destroyed during
# the course of the algorithm.
#
# Initially there are no non-trivial blossoms.
2024-06-25 21:01:52 +02:00
self.nontrivial_blossom: set[_NonTrivialBlossom] = set()
# "vertex_set_node[x]" represents the vertex "x" inside the
# union-find datastructure of its top-level blossom.
2024-06-29 22:06:28 +02:00
#
2024-06-30 21:38:12 +02:00
# Initially, each vertex belongs to its own trivial top-level blossom.
self.vertex_set_node = [b.vertex_set.insert(i, math.inf)
for (i, b) in enumerate(self.trivial_blossom)]
2024-05-26 20:55:54 +02:00
# All vertex duals are initialized to half the maximum edge weight.
#
2024-05-26 20:55:54 +02:00
# "start_vertex_dual_2x" is 2 times the initial vertex dual value.
#
# Pre-multiplication by 2 ensures that the values are integers
# if all edge weights are integers.
self.start_vertex_dual_2x = max(w for (_x, _y, w) in graph.edges)
# Every vertex has a variable in the dual LPP.
2024-05-25 23:18:15 +02:00
#
2024-05-26 20:55:54 +02:00
# The value of the dual variable changes through delta steps,
# but these changes are implemented as lazy updates.
#
# The true dual value of an S-vertex is
# (vertex_dual_2x[x] - delta_sum_2x) / 2
#
# The true dual value of a T-vertex is
# (vertex_dual_2x[x] + delta_sum_2x + B(x).vertex_dual_offset) / 2
2024-05-25 23:18:15 +02:00
#
2024-05-26 20:55:54 +02:00
# The true dual value of an unlabeled vertex is
# (vertex_dual_2x[x] + B(x).vertex_dual_offset) / 2
2024-05-25 23:18:15 +02:00
#
2024-05-26 20:55:54 +02:00
# Note that "vertex_dual_2x" is invariant under delta steps.
self.vertex_dual_2x: list[float]
self.vertex_dual_2x = num_vertex * [self.start_vertex_dual_2x]
# Running sum of applied delta steps times 2.
self.delta_sum_2x: float = 0
# Queue containing unlabeled top-level blossoms that have an edge to
2024-06-29 22:06:28 +02:00
# an S-blossom. The priority of a blossom is 2 times its least slack
# to an S blossom, plus 2 times the running sum of delta steps.
self.delta2_queue: PriorityQueue[_Blossom] = PriorityQueue()
# Queue containing edges between S-vertices in different top-level
# blossoms. The priority of an edge is its slack plus 2 times the
2024-05-26 00:07:16 +02:00
# running sum of delta steps.
self.delta3_queue: PriorityQueue[int] = PriorityQueue()
self.delta3_node: list[Optional[PriorityQueue.Node]]
self.delta3_node = [None for _e in graph.edges]
# Queue containing top-level non-trivial T-blossoms.
# The priority of a blossom is its dual plus 2 times the running
2024-05-26 00:07:16 +02:00
# sum of delta steps.
self.delta4_queue: PriorityQueue[_NonTrivialBlossom] = PriorityQueue()
# For each T-vertex or unlabeled vertex "x",
# "vertex_sedge_queue[x]" is a queue of edges between "x" and any
# S-vertex. The priority of an edge is 2 times its pseudo-slack.
2024-07-06 16:13:34 +02:00
self.vertex_sedge_queue: list[PriorityQueue[int]]
self.vertex_sedge_queue = [PriorityQueue() for _x in range(num_vertex)]
self.vertex_sedge_node: list[Optional[PriorityQueue.Node]]
self.vertex_sedge_node = [None for _e in graph.edges]
2023-04-06 20:50:06 +02:00
# Queue of S-vertices to be scanned.
self.scan_queue: list[int] = []
2024-05-25 23:18:15 +02:00
def __del__(self) -> None:
"""Delete reference cycles during cleanup of the matching context."""
2024-06-25 21:01:52 +02:00
for blossom in itertools.chain(self.trivial_blossom,
self.nontrivial_blossom):
2024-07-06 16:13:34 +02:00
blossom.parent = None
blossom.vertex_set.clear()
del blossom.vertex_set
2024-06-29 22:06:28 +02:00
#
# Least-slack edge tracking:
#
def edge_pseudo_slack_2x(self, e: int) -> float:
"""Return 2 times the pseudo-slack of the specified edge.
The pseudo-slack of an edge is related to its true slack, but
2024-06-29 22:06:28 +02:00
adjusted in a way that makes it invariant under delta steps.
If the edge connects two S-vertices in different top-level blossoms,
the true slack is the pseudo-slack minus 2 times the running sum
of delta steps.
If the edge connects an S-vertex to an unlabeled vertex,
the true slack is the pseudo-slack minus the running sum of delta
steps, plus the pending offset of the top-level blossom that contains
the unlabeled vertex.
"""
(x, y, w) = self.graph.edges[e]
return self.vertex_dual_2x[x] + self.vertex_dual_2x[y] - 2 * w
2024-06-29 22:06:28 +02:00
def delta2_add_edge(self, e: int, y: int, by: _Blossom) -> None:
"""Add edge "e" for delta2 tracking.
2024-06-29 22:06:28 +02:00
Edge "e" connects an S-vertex to a T-vertex or unlabeled vertex "y".
This function takes time O(log(n)).
"""
2024-06-29 22:06:28 +02:00
prio = self.edge_pseudo_slack_2x(e)
improved = (self.vertex_sedge_queue[y].empty()
or (self.vertex_sedge_queue[y].find_min().prio > prio))
2024-06-29 22:06:28 +02:00
# Insert edge in the S-edge queue of vertex "y".
assert self.vertex_sedge_node[e] is None
self.vertex_sedge_node[e] = self.vertex_sedge_queue[y].insert(prio, e)
2024-06-29 22:06:28 +02:00
# Continue if the new edge becomes the least-slack S-edge for "y".
if not improved:
return
2024-06-29 22:06:28 +02:00
# Update the priority of "y" in its UnionFindQueue.
self.vertex_set_node[y].set_prio(prio)
2024-06-29 22:06:28 +02:00
# If the blossom is unlabeled and the new edge becomes its least-slack
# S-edge, insert or update the blossom in the global delta2 queue.
2024-07-06 20:21:33 +02:00
if by.label == _LABEL_NONE:
prio += by.vertex_dual_offset
if by.delta2_node is None:
by.delta2_node = self.delta2_queue.insert(prio, by)
elif prio < by.delta2_node.prio:
self.delta2_queue.decrease_prio(by.delta2_node, prio)
2024-06-29 22:06:28 +02:00
def delta2_remove_edge(self, e: int, y: int, by: _Blossom) -> None:
"""Remove edge "e" from delta2 tracking.
This function is called if an S-vertex becomes unlabeled,
and edge "e" connects that vertex to vertex "y" which is a T-vertex
or unlabeled vertex.
This function takes time O(log(n)).
"""
vertex_sedge_node = self.vertex_sedge_node[e]
if vertex_sedge_node is not None:
# Delete edge from the S-edge queue of vertex "y".
vertex_sedge_queue = self.vertex_sedge_queue[y]
vertex_sedge_queue.delete(vertex_sedge_node)
self.vertex_sedge_node[e] = None
if vertex_sedge_queue.empty():
prio = math.inf
else:
prio = vertex_sedge_queue.find_min().prio
# If necessary, update the priority of "y" in its UnionFindQueue.
if prio > self.vertex_set_node[y].prio:
self.vertex_set_node[y].set_prio(prio)
if by.label == _LABEL_NONE:
# Update or delete the blossom in the global delta2 queue.
assert by.delta2_node is not None
prio = by.vertex_set.min_prio()
if prio < math.inf:
prio += by.vertex_dual_offset
if prio > by.delta2_node.prio:
self.delta2_queue.increase_prio(
by.delta2_node, prio)
else:
self.delta2_queue.delete(by.delta2_node)
by.delta2_node = None
def delta2_enable_blossom(self, blossom: _Blossom) -> None:
"""Enable delta2 tracking for "blossom".
This function is called when a blossom becomes an unlabeled top-level
blossom. If the blossom has at least one edge to an S-vertex,
the blossom will be inserted in the global delta2 queue.
This function takes time O(log(n)).
"""
assert blossom.delta2_node is None
prio = blossom.vertex_set.min_prio()
if prio < math.inf:
prio += blossom.vertex_dual_offset
blossom.delta2_node = self.delta2_queue.insert(prio, blossom)
def delta2_disable_blossom(self, blossom: _Blossom) -> None:
"""Disable delta2 tracking for "blossom".
The blossom will be removed from the global delta2 queue.
This function is called when a blossom stops being an unlabeled
top-level blossom.
This function takes time O(log(n)).
"""
if blossom.delta2_node is not None:
self.delta2_queue.delete(blossom.delta2_node)
blossom.delta2_node = None
def delta2_clear_vertex(self, x: int) -> None:
"""Clear delta2 tracking for vertex "x".
This function is called when "x" becomes an S-vertex.
It is assumed that the blossom containing "x" has already been
disabled for delta2 tracking.
2024-06-30 21:38:12 +02:00
This function takes time O(k + log(n)),
2024-06-29 22:06:28 +02:00
where "k" is the number of edges incident on "x".
"""
self.vertex_sedge_queue[x].clear()
for e in self.graph.adjacent_edges[x]:
self.vertex_sedge_node[e] = None
self.vertex_set_node[x].set_prio(math.inf)
def delta2_get_min_edge(self) -> tuple[int, float]:
"""Find the least-slack edge between any S-vertex and any unlabeled
vertex.
This function takes time O(log(n)).
Returns:
2024-06-29 22:06:28 +02:00
Tuple (edge_index, slack_2x) if there is an S-to-unlabeled edge,
or (-1, Inf) if there is no such edge.
"""
if self.delta2_queue.empty():
2024-06-29 22:06:28 +02:00
return (-1, math.inf)
delta2_node = self.delta2_queue.find_min()
blossom = delta2_node.data
prio = delta2_node.prio
slack_2x = prio - self.delta_sum_2x
assert blossom.parent is None
assert blossom.label == _LABEL_NONE
x = blossom.vertex_set.min_elem()
e = self.vertex_sedge_queue[x].find_min().data
return (e, slack_2x)
2024-06-29 22:06:28 +02:00
def delta3_add_edge(self, e: int) -> None:
"""Add edge "e" for delta3 tracking.
This function is called if a vertex becomes an S-vertex and edge "e"
connects it to an S-vertex in a different top-level blossom.
This function takes time O(log(n)).
"""
# The edge may already be in the delta3 queue, if it was previously
# discovered in the opposite direction.
if self.delta3_node[e] is None:
# Priority is edge slack plus 2 times the running sum of
# delta steps.
prio_2x = self.edge_pseudo_slack_2x(e)
if self.graph.integer_weights:
# If all edge weights are integers, the slack of
# any edge between S-vertices is also an integer.
assert prio_2x % 2 == 0
prio = prio_2x // 2
else:
prio = prio_2x / 2
self.delta3_node[e] = self.delta3_queue.insert(prio, e)
def delta3_remove_edge(self, e: int) -> None:
"""Remove edge "e" from delta3 tracking.
This function is called if a former S-vertex becomes unlabeled,
and edge "e" connects it to another S-vertex.
This function takes time O(log(n)).
"""
delta3_node = self.delta3_node[e]
if delta3_node is not None:
self.delta3_queue.delete(delta3_node)
self.delta3_node[e] = None
def delta3_get_min_edge(self) -> tuple[int, float]:
"""Find the least-slack edge between any pair of S-vertices in
different top-level blossoms.
This function takes time O(1 + k * log(n)),
where "k" is the number of intra-blossom edges removed from the queue.
Returns:
Tuple (edge_index, slack) if there is an S-to-S edge,
or (-1, Inf) if there is no suitable edge.
"""
while not self.delta3_queue.empty():
delta3_node = self.delta3_queue.find_min()
e = delta3_node.data
(x, y, _w) = self.graph.edges[e]
bx = self.vertex_set_node[x].find()
by = self.vertex_set_node[y].find()
assert (bx.label == _LABEL_S) and (by.label == _LABEL_S)
if bx is not by:
slack = delta3_node.prio - self.delta_sum_2x
return (e, slack)
# Reject edges between vertices within the same top-level blossom.
# Although intra-blossom edges are never inserted into the queue,
# existing edges in the queue may become intra-blossom when
# a new blossom is formed.
self.delta3_queue.delete(delta3_node)
self.delta3_node[e] = None
# If the queue is empty, no suitable edge exists.
return (-1, math.inf)
#
2024-06-30 09:11:58 +02:00
# Managing blossom labels:
#
2023-02-04 21:06:34 +01:00
2024-05-25 23:18:15 +02:00
def assign_blossom_label_s(self, blossom: _Blossom) -> None:
2024-06-30 21:38:12 +02:00
"""Change an unlabeled top-level blossom into an S-blossom.
For a blossom with "j" vertices and "k" incident edges,
this function takes time O(j * log(n) + k).
This function is called at most once per blossom per stage.
It therefore takes total time O(n * log(n) + m) per stage.
"""
2024-06-30 09:11:58 +02:00
2024-05-26 00:07:16 +02:00
assert blossom.parent is None
assert blossom.label == _LABEL_NONE
blossom.label = _LABEL_S
2024-05-25 23:18:15 +02:00
2024-06-30 09:11:58 +02:00
# Labeled blossoms must not be in the delta2 queue.
2024-06-29 22:06:28 +02:00
self.delta2_disable_blossom(blossom)
2024-06-30 09:11:58 +02:00
# Adjust for lazy updating of S-blossom dual variables.
#
# The true dual value of an unlabeled top-level blossom is
# blossom.dual_var
#
# while the true dual value of a top-level S-blossom is
# blossom.dual_var + ctx.delta_sum_2x
#
# The value of blossom.dual_var must be adjusted accordingly
# when the blossom changes from unlabeled to S-blossom.
#
2024-05-26 00:07:16 +02:00
if isinstance(blossom, _NonTrivialBlossom):
blossom.dual_var -= self.delta_sum_2x
2024-05-25 23:18:15 +02:00
2024-06-30 09:11:58 +02:00
# Apply pending updates to vertex dual variables and prepare
# for lazy updating of S-vertex dual variables.
#
# For S-blossoms, blossom.vertex_dual_offset is always 0.
#
# Furthermore, the true dual value of an unlabeled vertex is
# (vertex_dual_2x[x] + blossom.vertex_dual_offset) / 2
#
# while the true dual value of an S-vertex is
# (vertex_dual_2x[x] - delta_sum_2x) / 2
#
# The value of vertex_dual_2x must be adjusted accordingly
# when vertices change from unlabeled to S-vertex.
#
vertex_dual_fixup = self.delta_sum_2x + blossom.vertex_dual_offset
blossom.vertex_dual_offset = 0
2024-06-30 09:11:58 +02:00
vertices = blossom.vertices()
2024-05-26 20:55:54 +02:00
for x in vertices:
self.vertex_dual_2x[x] += vertex_dual_fixup
2024-06-30 09:11:58 +02:00
# S-vertices do not keep track of potential delta2 edges.
2024-06-29 22:06:28 +02:00
self.delta2_clear_vertex(x)
2024-06-30 09:11:58 +02:00
# Add the new S-vertices to the scan queue.
self.scan_queue.extend(vertices)
2024-05-25 23:18:15 +02:00
def assign_blossom_label_t(self, blossom: _Blossom) -> None:
2024-06-30 21:38:12 +02:00
"""Change an unlabeled top-level blossom into a T-blossom.
This function takes time O(log(n)).
"""
2024-05-25 23:18:15 +02:00
assert blossom.parent is None
assert blossom.label == _LABEL_NONE
blossom.label = _LABEL_T
2024-06-30 09:11:58 +02:00
# Labeled blossoms must not be in the delta2 queue.
2024-06-29 22:06:28 +02:00
self.delta2_disable_blossom(blossom)
2024-05-25 23:18:15 +02:00
if isinstance(blossom, _NonTrivialBlossom):
2024-05-26 11:32:02 +02:00
2024-06-30 09:11:58 +02:00
# Adjust for lazy updating of T-blossom dual variables.
#
# The true dual value of an unlabeled top-level blossom is
# blossom.dual_var
#
# while the true dual value of a top-level T-blossom is
# blossom.dual_var - ctx.delta_sum_2x
#
# The value of blossom.dual_var must be adjusted accordingly
# when the blossom changes from unlabeled to S-blossom.
#
2024-05-26 11:32:02 +02:00
blossom.dual_var += self.delta_sum_2x
2024-06-30 09:11:58 +02:00
# Top-level T-blossoms are tracked in the delta4 queue.
2024-05-25 23:18:15 +02:00
assert blossom.delta4_node is None
2024-05-26 11:32:02 +02:00
blossom.delta4_node = self.delta4_queue.insert(blossom.dual_var,
blossom)
2024-05-25 23:18:15 +02:00
# Prepare for lazy updating of T-vertex dual variables.
2024-06-30 09:11:58 +02:00
#
# The true dual value of an unlabeled vertex is
# (vertex_dual_2x[x] + blossom.vertex_dual_offset) / 2
#
# while the true dual value of a T-vertex is
# (vertex_dual_2x[x] + blossom.vertex_dual_offset + delta_sum_2x) / 2
#
# The value of blossom.vertex_dual_offset must be adjusted accordingly
# when the blossom changes from unlabeled to T-blossom.
#
blossom.vertex_dual_offset -= self.delta_sum_2x
2024-05-26 20:55:54 +02:00
2024-06-30 09:11:58 +02:00
def remove_blossom_label_s(self, blossom: _Blossom) -> None:
"""Change a top-level S-blossom into an unlabeled blossom.
For a blossom with "j" vertices and "k" incident edges,
this function takes time O((j + k) * log(n)).
This function is called at most once per blossom per stage.
It therefore takes total time O((n + m) * log(n)) per stage.
"""
assert blossom.parent is None
assert blossom.label == _LABEL_S
blossom.label = _LABEL_NONE
# Unwind lazy delta updates to the S-blossom dual variable.
if isinstance(blossom, _NonTrivialBlossom):
blossom.dual_var += self.delta_sum_2x
assert blossom.vertex_dual_offset == 0
vertex_dual_fixup = -self.delta_sum_2x
edges = self.graph.edges
adjacent_edges = self.graph.adjacent_edges
for x in blossom.vertices():
# Unwind lazy delta updates to S-vertex dual variables.
self.vertex_dual_2x[x] += vertex_dual_fixup
# Scan the incident edges of all vertices in the blossom.
for e in adjacent_edges[x]:
(p, q, _w) = edges[e]
y = p if p != x else q
# If this edge is in the delta3 queue, remove it.
# Only edges between S-vertices are tracked for delta3,
# and vertex "x" is no longer an S-vertex.
self.delta3_remove_edge(e)
by = self.vertex_set_node[y].find()
if by.label == _LABEL_S:
# Edge "e" connects unlabeled vertex "x" to S-vertex "y".
# It must be tracked for delta2 via vertex "x".
self.delta2_add_edge(e, x, blossom)
else:
# Edge "e" connects former S-vertex "x" to T-vertex
# or unlabeled vertex "y". That implies this edge was
# tracked for delta2 via vertex "y", but it must be
# removed now.
self.delta2_remove_edge(e, y, by)
2024-05-25 23:18:15 +02:00
def remove_blossom_label_t(self, blossom: _Blossom) -> None:
2024-06-30 21:38:12 +02:00
"""Change a top-level T-blossom into an unlabeled blossom.
This function takes time O(log(n)).
"""
2024-05-25 23:18:15 +02:00
assert blossom.parent is None
assert blossom.label == _LABEL_T
blossom.label = _LABEL_NONE
if isinstance(blossom, _NonTrivialBlossom):
2024-05-26 11:32:02 +02:00
2024-06-30 09:11:58 +02:00
# Unlabeled blossoms are not tracked in the delta4 queue.
2024-05-25 23:18:15 +02:00
assert blossom.delta4_node is not None
self.delta4_queue.delete(blossom.delta4_node)
blossom.delta4_node = None
2024-06-30 09:11:58 +02:00
# Unwind lazy updates to the T-blossom dual variable.
2024-05-26 11:32:02 +02:00
blossom.dual_var -= self.delta_sum_2x
# Unwind lazy updates of T-vertex dual variables.
blossom.vertex_dual_offset += self.delta_sum_2x
2024-05-26 20:55:54 +02:00
2024-06-30 09:11:58 +02:00
# Enable unlabeled top-level blossom for delta2 tracking.
self.delta2_enable_blossom(blossom)
2024-06-30 09:11:58 +02:00
def change_s_blossom_to_subblossom(self, blossom: _Blossom) -> None:
2024-06-30 21:38:12 +02:00
"""Change a top-level S-blossom into an S-subblossom.
This function takes time O(1).
"""
2024-06-30 09:11:58 +02:00
assert blossom.parent is None
assert blossom.label == _LABEL_S
blossom.label = _LABEL_NONE
2024-06-30 09:11:58 +02:00
# Unwind lazy delta updates to the S-blossom dual variable.
if isinstance(blossom, _NonTrivialBlossom):
blossom.dual_var += self.delta_sum_2x
2024-06-30 09:11:58 +02:00
#
# General support routines:
#
2024-05-26 20:55:54 +02:00
def reset_blossom_label(self, blossom: _Blossom) -> None:
"""Remove blossom label."""
2024-05-26 20:55:54 +02:00
assert blossom.parent is None
assert blossom.label != _LABEL_NONE
2024-05-26 20:55:54 +02:00
if blossom.label == _LABEL_S:
2024-06-30 09:11:58 +02:00
self.remove_blossom_label_s(blossom)
2024-07-06 20:21:33 +02:00
else:
self.remove_blossom_label_t(blossom)
def _check_alternating_tree_consistency(self) -> None:
"""TODO -- remove this function, only for debugging"""
2024-06-25 21:01:52 +02:00
for blossom in itertools.chain(self.trivial_blossom,
self.nontrivial_blossom):
if (blossom.parent is None) and (blossom.label != _LABEL_NONE):
assert blossom.tree_blossoms is not None
assert blossom in blossom.tree_blossoms
if blossom.tree_edge is not None:
bx = self.vertex_set_node[blossom.tree_edge[0]].find()
by = self.vertex_set_node[blossom.tree_edge[1]].find()
assert bx.tree_blossoms is blossom.tree_blossoms
assert by.tree_blossoms is blossom.tree_blossoms
else:
assert blossom.tree_edge is None
assert blossom.tree_blossoms is None
def remove_alternating_tree(
self,
tree_blossoms: set[_Blossom]
) -> None:
"""Reset the alternating tree consisting of the specified blossoms.
Marks the blossoms as unlabeled.
Updates delta tracking accordingly.
This function takes time O((n + m) * log(n)).
"""
for blossom in tree_blossoms:
assert blossom.label != _LABEL_NONE
assert blossom.tree_blossoms is tree_blossoms
self.reset_blossom_label(blossom)
blossom.tree_edge = None
blossom.tree_blossoms = None
2023-02-07 11:23:42 +01:00
def trace_alternating_paths(self, x: int, y: int) -> _AlternatingPath:
"""Trace back through the alternating trees from vertices "x" and "y".
2023-02-04 21:06:34 +01:00
2023-02-06 22:33:26 +01:00
If both vertices are part of the same alternating tree, this function
discovers a new blossom. In this case it returns an alternating path
through the blossom that starts and ends in the same sub-blossom.
2023-02-04 21:06:34 +01:00
2023-02-06 22:33:26 +01:00
If the vertices are part of different alternating trees, this function
discovers an augmenting path. In this case it returns an alternating
path that starts and ends in an unmatched vertex.
2023-02-04 21:06:34 +01:00
This function takes time O(k * log(n)) to discover a blossom,
where "k" is the number of sub-blossoms,
or time O(n * log(n)) to discover an augmenting path.
2023-02-04 21:06:34 +01:00
2023-02-06 22:33:26 +01:00
Returns:
Alternating path as an ordered list of edges between top-level
blossoms.
"""
2023-02-04 21:06:34 +01:00
marked_blossoms: list[_Blossom] = []
2023-02-06 22:33:26 +01:00
2023-02-07 11:23:42 +01:00
# "xedges" is a list of edges used while tracing from "x".
# "yedges" is a list of edges used while tracing from "y".
2023-02-08 19:21:34 +01:00
# Pre-load the edge (x, y) on both lists.
xedges: list[tuple[int, int]] = [(x, y)]
yedges: list[tuple[int, int]] = [(y, x)]
2023-02-06 22:33:26 +01:00
2023-02-08 19:21:34 +01:00
# "first_common" is the first common ancestor of "x" and "y"
# in the alternating tree, or None if there is no common ancestor.
first_common: Optional[_Blossom] = None
2023-02-06 22:33:26 +01:00
2023-02-07 11:23:42 +01:00
# Alternate between tracing the path from "x" and the path from "y".
2023-02-06 22:33:26 +01:00
# This ensures that the search time is bounded by the size of the
# newly found blossom.
2023-02-07 11:23:42 +01:00
while x != -1 or y != -1:
2023-02-06 22:33:26 +01:00
# Check if we found a common ancestor.
bx = self.vertex_set_node[x].find()
if bx.marker:
2023-02-07 11:23:42 +01:00
first_common = bx
2023-02-06 22:33:26 +01:00
break
2023-02-04 21:06:34 +01:00
2023-02-06 22:33:26 +01:00
# Mark blossom as a potential common ancestor.
bx.marker = True
2023-02-07 11:23:42 +01:00
marked_blossoms.append(bx)
2023-02-04 21:06:34 +01:00
2023-02-06 22:33:26 +01:00
# Track back through the link in the alternating tree.
if bx.tree_edge is None:
2023-02-06 22:33:26 +01:00
# Reached the root of this alternating tree.
2023-02-07 11:23:42 +01:00
x = -1
else:
xedges.append(bx.tree_edge)
x = bx.tree_edge[0]
2023-02-04 21:06:34 +01:00
2023-02-07 11:23:42 +01:00
# Swap "x" and "y" to alternate between paths.
if y != -1:
(x, y) = (y, x)
(xedges, yedges) = (yedges, xedges)
2023-02-04 21:06:34 +01:00
2023-02-06 22:33:26 +01:00
# Remove all markers we placed.
for b in marked_blossoms:
b.marker = False
2023-02-04 21:06:34 +01:00
2023-02-06 22:33:26 +01:00
# If we found a common ancestor, trim the paths so they end there.
if first_common is not None:
assert self.vertex_set_node[xedges[-1][0]].find() is first_common
while (self.vertex_set_node[yedges[-1][0]].find()
is not first_common):
2023-02-07 11:23:42 +01:00
yedges.pop()
2023-02-04 21:06:34 +01:00
2023-02-06 22:33:26 +01:00
# Fuse the two paths.
2023-02-08 19:21:34 +01:00
# Flip the order of one path, and flip the edge tuples in the other
# path to obtain a continuous path with correctly ordered edge tuples.
# Skip the duplicate edge in one of the paths.
path_edges = xedges[::-1] + [(y, x) for (x, y) in yedges[1:]]
2023-02-04 21:06:34 +01:00
2023-02-06 22:33:26 +01:00
# Any S-to-S alternating path must have odd length.
assert len(path_edges) % 2 == 1
2023-02-04 21:06:34 +01:00
return _AlternatingPath(edges=path_edges,
is_cycle=(first_common is not None))
2023-02-05 19:33:13 +01:00
#
# Merge and expand blossoms:
#
2023-02-06 22:33:26 +01:00
def make_blossom(self, path: _AlternatingPath) -> None:
"""Create a new blossom from an alternating cycle.
2023-02-05 19:33:13 +01:00
2023-02-06 22:33:26 +01:00
Assign label S to the new blossom.
Relabel all T-sub-blossoms as S and add their vertices to the queue.
2023-02-05 19:33:13 +01:00
2024-06-30 21:38:12 +02:00
This function takes total time O((n + m) * log(n)) per stage.
2023-02-06 22:33:26 +01:00
"""
2023-02-04 21:06:34 +01:00
2023-02-06 22:33:26 +01:00
# Check that the path is odd-length.
assert len(path.edges) % 2 == 1
assert len(path.edges) >= 3
2023-02-04 21:06:34 +01:00
2023-02-06 22:33:26 +01:00
# Construct the list of sub-blossoms (current top-level blossoms).
subblossoms = [self.vertex_set_node[x].find() for (x, y) in path.edges]
2023-02-04 21:06:34 +01:00
2023-02-06 22:33:26 +01:00
# Check that the path is cyclic.
2024-05-25 23:18:15 +02:00
# Note the path will not always start and end with the same _vertex_,
2023-02-06 22:33:26 +01:00
# but it must start and end in the same _blossom_.
subblossoms_next = [self.vertex_set_node[y].find()
for (x, y) in path.edges]
2023-02-06 22:33:26 +01:00
assert subblossoms[0] == subblossoms_next[-1]
assert subblossoms[1:] == subblossoms_next[:-1]
2023-02-04 21:06:34 +01:00
2024-05-25 23:18:15 +02:00
# Blossom must start and end with an S-sub-blossom.
assert subblossoms[0].label == _LABEL_S
2024-05-26 00:07:16 +02:00
# Remove blossom labels.
# Mark vertices inside former T-blossoms as S-vertices.
2024-05-25 23:18:15 +02:00
for sub in subblossoms:
2024-06-30 09:11:58 +02:00
if sub.label == _LABEL_T:
2024-05-25 23:18:15 +02:00
self.remove_blossom_label_t(sub)
2024-06-30 09:11:58 +02:00
self.assign_blossom_label_s(sub)
self.change_s_blossom_to_subblossom(sub)
2024-05-25 23:18:15 +02:00
2023-02-06 22:33:26 +01:00
# Create the new blossom object.
blossom = _NonTrivialBlossom(subblossoms, path.edges)
2023-02-06 22:33:26 +01:00
2024-06-30 09:11:58 +02:00
# Assign label S to the new blossom.
blossom.label = _LABEL_S
2024-06-30 09:11:58 +02:00
# Prepare for lazy updating of S-blossom dual variable.
blossom.dual_var = -self.delta_sum_2x
# Link the new blossom to the alternating tree.
tree_blossoms = subblossoms[0].tree_blossoms
assert tree_blossoms is not None
2024-05-25 23:18:15 +02:00
blossom.tree_edge = subblossoms[0].tree_edge
blossom.tree_blossoms = tree_blossoms
tree_blossoms.add(blossom)
2024-05-25 23:18:15 +02:00
2024-06-30 09:11:58 +02:00
# Add to the list of blossoms.
2024-06-25 21:01:52 +02:00
self.nontrivial_blossom.add(blossom)
2023-02-06 22:33:26 +01:00
# Link the subblossoms to the their new parent.
for sub in subblossoms:
sub.parent = blossom
2023-02-06 22:33:26 +01:00
# Remove subblossom from the alternating tree.
sub.tree_edge = None
sub.tree_blossoms = None
tree_blossoms.remove(sub)
# Merge union-find structures.
blossom.vertex_set.merge([sub.vertex_set for sub in subblossoms])
2023-02-08 19:24:14 +01:00
@staticmethod
2023-02-06 22:33:26 +01:00
def find_path_through_blossom(
blossom: _NonTrivialBlossom,
sub: _Blossom
) -> tuple[list[_Blossom], list[tuple[int, int]]]:
"""Construct a path with an even number of edges through the
specified blossom, from sub-blossom "sub" to the base of "blossom".
2023-02-06 22:33:26 +01:00
Return:
Tuple (nodes, edges).
"""
2023-02-04 21:06:34 +01:00
2023-02-06 22:33:26 +01:00
# Walk around the blossom from "sub" to its base.
p = blossom.subblossoms.index(sub)
2023-05-12 18:12:25 +02:00
if p % 2 == 0:
# Walk backwards around the blossom.
# Flip edges from (i,j) to (j,i) to make them fit
# in the path from "sub" to base.
nodes = blossom.subblossoms[p::-1]
edges = [(j, i) for (i, j) in blossom.edges[:p][::-1]]
else:
# Walk forward around the blossom.
nodes = blossom.subblossoms[p:] + blossom.subblossoms[0:1]
edges = blossom.edges[p:]
2023-02-06 22:33:26 +01:00
assert len(edges) % 2 == 0
assert len(nodes) % 2 == 1
2023-02-06 22:33:26 +01:00
return (nodes, edges)
2024-06-30 09:11:58 +02:00
def expand_unlabeled_blossom(self, blossom: _NonTrivialBlossom) -> None:
"""Expand the specified unlabeled blossom.
2023-02-06 22:33:26 +01:00
2024-06-30 21:38:12 +02:00
This function takes total time O(n * log(n)) per stage.
2023-02-06 22:33:26 +01:00
"""
2023-02-04 21:06:34 +01:00
assert blossom.parent is None
2024-06-30 09:11:58 +02:00
assert blossom.label == _LABEL_NONE
2024-06-30 09:11:58 +02:00
# Remove blossom from the delta2 queue.
self.delta2_disable_blossom(blossom)
# Split union-find structure.
blossom.vertex_set.split()
2024-05-26 20:55:54 +02:00
# Prepare to push lazy delta updates down to the sub-blossoms.
2024-06-30 09:11:58 +02:00
vertex_dual_offset = blossom.vertex_dual_offset
2024-05-26 20:55:54 +02:00
blossom.vertex_dual_offset = 0
2023-02-06 22:33:26 +01:00
# Convert sub-blossoms into top-level blossoms.
for sub in blossom.subblossoms:
assert sub.label == _LABEL_NONE
sub.parent = None
assert sub.vertex_dual_offset == 0
2024-06-30 09:11:58 +02:00
sub.vertex_dual_offset = vertex_dual_offset
2024-06-29 22:06:28 +02:00
self.delta2_enable_blossom(sub)
2023-02-06 22:33:26 +01:00
2024-07-06 16:13:34 +02:00
# Avoid leaking a reference cycle.
del blossom.vertex_set
2024-06-30 09:11:58 +02:00
# Delete the expanded blossom.
self.nontrivial_blossom.remove(blossom)
def expand_t_blossom(self, blossom: _NonTrivialBlossom) -> None:
"""Expand the specified T-blossom.
2024-06-30 21:38:12 +02:00
This function takes total time O(n * log(n) + m) per stage.
2024-06-30 09:11:58 +02:00
"""
assert blossom.parent is None
assert blossom.label == _LABEL_T
assert blossom.delta2_node is None
# Remove blossom from its alternating tree.
tree_blossoms = blossom.tree_blossoms
assert tree_blossoms is not None
tree_blossoms.remove(blossom)
# Remove label T.
self.remove_blossom_label_t(blossom)
# Expand the now-unlabeled blossom.
self.expand_unlabeled_blossom(blossom)
# The expanded blossom was part of an alternating tree, linked to
2023-02-06 22:33:26 +01:00
# a parent node in the tree via one of its subblossoms, and linked to
# a child node of the tree via the base vertex.
# We must reconstruct this part of the alternating tree, which will
# now run via sub-blossoms of the expanded blossom.
# Find the sub-blossom that is attached to the parent node in
# the alternating tree.
2023-02-10 12:58:47 +01:00
assert blossom.tree_edge is not None
(x, y) = blossom.tree_edge
sub = self.vertex_set_node[y].find()
2023-02-04 21:06:34 +01:00
2023-02-06 22:33:26 +01:00
# Assign label T to that sub-blossom.
2024-05-25 23:18:15 +02:00
self.assign_blossom_label_t(sub)
sub.tree_edge = blossom.tree_edge
sub.tree_blossoms = tree_blossoms
tree_blossoms.add(sub)
2023-02-04 21:06:34 +01:00
2023-02-06 22:33:26 +01:00
# 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) = self.find_path_through_blossom(blossom,
sub)
2023-02-04 21:06:34 +01:00
2023-02-06 22:33:26 +01:00
for p in range(0, len(path_edges), 2):
#
2023-02-07 11:23:42 +01:00
# (p) ==(y,x)== (p+1) ----- (p+2)
2023-02-06 22:33:26 +01:00
# T S T
#
# path_nodes[p] has already been labeled T.
# We now assign labels to path_nodes[p+1] and path_nodes[p+2].
2023-02-04 21:06:34 +01:00
2023-02-06 22:33:26 +01:00
# Assign label S to path_nodes[p+1].
2023-02-07 11:23:42 +01:00
(y, x) = path_edges[p]
self.extend_tree_t_to_s(x)
2023-02-04 21:06:34 +01:00
# Assign label T to path_nodes[i+2] and attach it
# to path_nodes[p+1].
2023-02-06 22:33:26 +01:00
sub = path_nodes[p+2]
2024-05-25 23:18:15 +02:00
self.assign_blossom_label_t(sub)
sub.tree_edge = path_edges[p+1]
sub.tree_blossoms = tree_blossoms
tree_blossoms.add(sub)
2023-02-04 21:06:34 +01:00
#
# Augmenting:
#
2023-02-04 21:06:34 +01:00
2023-02-07 16:15:59 +01:00
def augment_blossom_rec(
self,
blossom: _NonTrivialBlossom,
sub: _Blossom,
stack: list[tuple[_NonTrivialBlossom, _Blossom]]
2023-02-07 16:15:59 +01:00
) -> None:
"""Augment along an alternating path through the specified blossom,
2023-02-07 16:15:59 +01:00
from sub-blossom "sub" to the base vertex of the blossom.
2023-02-04 21:06:34 +01:00
2023-02-07 16:15:59 +01:00
Modify the blossom to reflect that sub-blossom "sub" contains
the base vertex after augmenting.
2023-02-04 21:06:34 +01:00
Mark any sub-blossoms on the alternating path for recursive
augmentation, except for sub-blossom "sub" which has already been
augmented. Use the stack instead of making direct recursive calls.
2023-02-07 16:15:59 +01:00
"""
2023-02-04 21:06:34 +01:00
# Walk through the blossom from "sub" to the base vertex.
(path_nodes, path_edges) = self.find_path_through_blossom(blossom,
sub)
2023-02-07 16:15:59 +01:00
for p in range(0, len(path_edges), 2):
# Before augmentation:
# path_nodes[p] is matched to path_nodes[p+1]
#
# (p) ===== (p+1) ---(x,y)--- (p+2)
2023-02-07 16:15:59 +01:00
#
# After augmentation:
# path_nodes[p+1] matched to path_nodes[p+2] via edge (i,j)
#
# (p) ----- (p+1) ===(x,y)=== (p+2)
2023-02-07 16:15:59 +01:00
#
# Pull the edge (x, y) into the matching.
2023-02-07 16:15:59 +01:00
(x, y) = path_edges[p+1]
self.vertex_mate[x] = y
self.vertex_mate[y] = x
# Augment through the subblossoms touching the edge (x, y).
2023-02-07 16:15:59 +01:00
# Nothing needs to be done for trivial subblossoms.
bx = path_nodes[p+1]
if isinstance(bx, _NonTrivialBlossom):
stack.append((bx, self.trivial_blossom[x]))
2023-02-04 21:06:34 +01:00
2023-02-07 16:15:59 +01:00
by = path_nodes[p+2]
if isinstance(by, _NonTrivialBlossom):
stack.append((by, self.trivial_blossom[y]))
2023-02-07 16:15:59 +01:00
# Rotate the subblossom list so the new base ends up in position 0.
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
# its augmentation has already finished.
blossom.base_vertex = sub.base_vertex
2023-02-04 21:06:34 +01:00
def augment_blossom(
self,
blossom: _NonTrivialBlossom,
sub: _Blossom
) -> None:
"""Augment along an alternating path through the specified blossom,
2023-02-07 11:23:42 +01:00
from sub-blossom "sub" to the base vertex of the blossom.
2023-02-04 21:06:34 +01:00
Recursively augment any sub-blossoms on the alternating path.
2023-02-06 22:33:26 +01:00
This function takes time O(n).
"""
# Use an explicit stack to avoid deep recursion.
stack = [(blossom, sub)]
2023-02-06 22:33:26 +01:00
while stack:
2023-02-07 16:15:59 +01:00
(outer_blossom, sub) = stack.pop()
assert sub.parent is not None
blossom = sub.parent
2023-02-06 22:33:26 +01:00
if blossom != outer_blossom:
2023-02-07 16:15:59 +01:00
# Sub-blossom "sub" is an indirect (nested) child of
# the "outer_blossom" we are supposed to be augmenting.
2023-02-06 22:33:26 +01:00
#
# "blossom" is the direct parent of "sub".
# Let's first augment "blossom" from "sub" to its base vertex.
# Then continue by augmenting the parent of "blossom",
# from "blossom" to its base vertex, and so on until we
# get to the "outer_blossom".
2023-02-06 22:33:26 +01:00
#
# Set up to continue augmenting through the parent of
# "blossom".
stack.append((outer_blossom, blossom))
2023-02-06 22:33:26 +01:00
# Augment "blossom" from "sub" to the base vertex.
self.augment_blossom_rec(blossom, sub, stack)
2023-02-04 21:06:34 +01:00
2023-02-06 22:33:26 +01:00
def augment_matching(self, path: _AlternatingPath) -> None:
"""Augment the matching through the specified augmenting path.
2023-02-04 21:06:34 +01:00
2023-02-06 22:33:26 +01:00
This function takes time O(n).
"""
2023-02-04 21:06:34 +01:00
2023-02-06 22:33:26 +01:00
# Check that the augmenting path starts and ends in
# an unmatched vertex or a blossom with unmatched base.
assert len(path.edges) % 2 == 1
2023-02-07 11:23:42 +01:00
for x in (path.edges[0][0], path.edges[-1][1]):
b = self.vertex_set_node[x].find()
assert self.vertex_mate[b.base_vertex] == -1
2023-02-04 21:06:34 +01:00
2023-02-06 22:33:26 +01:00
# The augmenting path looks like this:
#
# (unmatched) ---- (B) ==== (B) ---- (B) ==== (B) ---- (unmatched)
#
# The first and last vertex (or blossom) of the path are unmatched
# (or have unmatched base vertex). After augmenting, those vertices
# will be matched. All matched edges on the path become unmatched,
# and unmatched edges become matched.
#
# This loop walks along the edges of this path that were not matched
# before augmenting.
2023-02-07 11:23:42 +01:00
for (x, y) in path.edges[0::2]:
2023-02-04 21:06:34 +01:00
2023-02-06 22:33:26 +01:00
# Augment the non-trivial blossoms on either side of this edge.
# No action is necessary for trivial blossoms.
bx = self.vertex_set_node[x].find()
if isinstance(bx, _NonTrivialBlossom):
self.augment_blossom(bx, self.trivial_blossom[x])
2023-02-04 21:06:34 +01:00
by = self.vertex_set_node[y].find()
if isinstance(by, _NonTrivialBlossom):
self.augment_blossom(by, self.trivial_blossom[y])
2023-02-04 21:06:34 +01:00
2023-02-06 22:33:26 +01:00
# Pull the edge into the matching.
2023-02-07 11:23:42 +01:00
self.vertex_mate[x] = y
self.vertex_mate[y] = x
2023-02-04 21:06:34 +01:00
#
# Alternating tree:
#
def extend_tree_t_to_s(self, x: int) -> None:
2023-02-07 11:23:42 +01:00
"""Assign label S to the unlabeled blossom that contains vertex "x".
2023-02-04 21:06:34 +01:00
The newly labeled S-blossom is added to the alternating tree
via its matched edge. All vertices in the newly labeled S-blossom
are added to the scan queue.
2023-02-04 21:06:34 +01:00
Preconditions:
- "x" is a vertex in an unlabeled blossom.
- "x" is matched to a T-vertex via a tight edge.
2023-02-06 22:33:26 +01:00
"""
2023-02-04 21:06:34 +01:00
# Assign label S to the blossom that contains vertex "x".
bx = self.vertex_set_node[x].find()
2024-05-25 23:18:15 +02:00
self.assign_blossom_label_s(bx)
2023-02-04 21:06:34 +01:00
# Vertex "x" is matched to T-vertex "y".
2023-02-07 11:23:42 +01:00
y = self.vertex_mate[x]
assert y != -1
2023-02-04 21:06:34 +01:00
by = self.vertex_set_node[y].find()
assert by.label == _LABEL_T
assert by.tree_blossoms is not None
2023-02-05 19:33:13 +01:00
# Attach the blossom that contains "x" to the alternating tree.
bx.tree_edge = (y, x)
bx.tree_blossoms = by.tree_blossoms
bx.tree_blossoms.add(bx)
def extend_tree_s_to_t(self, x: int, y: int) -> None:
2023-02-07 11:23:42 +01:00
"""Assign label T to the unlabeled blossom that contains vertex "y".
2023-02-05 19:33:13 +01:00
The newly labeled T-blossom is added to the alternating tree.
Directly afterwards, label S is assigned to the blossom that has
a matched edge to the base of the newly labeled T-blossom, and
that newly labeled S-blossom is also added to the alternating tree.
2023-02-05 19:33:13 +01:00
2023-02-06 22:33:26 +01:00
Preconditions:
2023-02-07 11:23:42 +01:00
- "x" is an S-vertex.
- "y" is a vertex in an unlabeled blossom with a matched base vertex.
2023-02-07 11:23:42 +01:00
- There is a tight edge between vertices "x" and "y".
2023-02-06 22:33:26 +01:00
"""
2023-02-05 19:33:13 +01:00
bx = self.vertex_set_node[x].find()
by = self.vertex_set_node[y].find()
assert bx.label == _LABEL_S
# Expand zero-dual blossoms before assigning label T.
while isinstance(by, _NonTrivialBlossom) and (by.dual_var == 0):
self.expand_unlabeled_blossom(by)
by = self.vertex_set_node[y].find()
# Assign label T to the unlabeled blossom.
2024-05-25 23:18:15 +02:00
self.assign_blossom_label_t(by)
by.tree_edge = (x, y)
by.tree_blossoms = bx.tree_blossoms
assert by.tree_blossoms is not None
by.tree_blossoms.add(by)
2023-02-05 19:33:13 +01:00
2023-05-10 20:54:29 +02:00
# Assign label S to the blossom that is mated to the T-blossom.
z = self.vertex_mate[by.base_vertex]
2023-02-07 11:23:42 +01:00
assert z != -1
self.extend_tree_t_to_s(z)
2023-02-05 19:33:13 +01:00
def add_s_to_s_edge(self, x: int, y: int) -> bool:
2023-02-07 11:23:42 +01:00
"""Add the edge between S-vertices "x" and "y".
2023-02-05 19:33:13 +01:00
2023-02-06 22:33:26 +01:00
If the edge connects blossoms that are part of the same alternating
tree, this function creates a new S-blossom and returns False.
2023-02-05 19:33:13 +01:00
2023-02-06 22:33:26 +01:00
If the edge connects two different alternating trees, an augmenting
path has been discovered. This function then augments the matching
and returns True. Labels are removed from blossoms that belonged
to the two alternating trees involved in the matching. All other
alternating trees and labels are preserved.
Preconditions:
- "x" and "y" are S-vertices in different top-level blossoms.
- There is a tight edge between vertices "x" and "y".
2023-02-05 19:33:13 +01:00
2023-02-06 22:33:26 +01:00
Returns:
True if the matching was augmented; otherwise False.
2023-02-06 22:33:26 +01:00
"""
2023-02-05 19:33:13 +01:00
bx = self.vertex_set_node[x].find()
by = self.vertex_set_node[y].find()
assert bx.label == _LABEL_S
assert by.label == _LABEL_S
assert bx is not by
2023-02-07 11:23:42 +01:00
# Trace back through the alternating trees from "x" and "y".
path = self.trace_alternating_paths(x, y)
2023-02-06 22:33:26 +01:00
assert bx.tree_blossoms is not None
assert by.tree_blossoms is not None
if bx.tree_blossoms is by.tree_blossoms:
# Both blossoms belong to the same alternating tree.
# This implies that the alternating path is a cycle.
# The path will be used to create a new blossom.
assert path.is_cycle
2023-02-06 22:33:26 +01:00
self.make_blossom(path)
return False
2023-02-05 19:33:13 +01:00
else:
# The blossoms belong to different alternating trees.
# This implies that the alternating path is an augmenting
# path between two unlabeled vertices.
# The path will be used to augment the matching.
# Delete the two alternating trees on the augmenting path.
# The blossoms in those trees become unlabeled.
self.remove_alternating_tree(bx.tree_blossoms)
self.remove_alternating_tree(by.tree_blossoms)
# Augment the matching.
assert not path.is_cycle
self.augment_matching(path)
return True
2023-02-05 19:33:13 +01:00
def scan_new_s_vertices(self) -> None:
"""Scan the incident edges of newly labeled S-vertices.
Edges are added to delta2 tracking or delta3 tracking depending
on the state of the vertex on the opposite side of the edge.
2023-02-04 21:06:34 +01:00
This function does not yet use the edges to extend the alternating
tree or find blossoms or augmenting paths, even if the edges
are tight. If there are such tight edges, they will be used later
through zero-delta steps.
If there are "j" new S-vertices with a total of "k" incident edges,
this function takes time O((j + k) * log(n)).
Since each vertex can become an S-vertex at most once per stage,
this function takes total time O((n + m) * log(n)) per stage.
2023-02-06 22:33:26 +01:00
"""
2023-02-04 21:06:34 +01:00
2023-02-06 22:33:26 +01:00
edges = self.graph.edges
adjacent_edges = self.graph.adjacent_edges
2023-02-04 21:06:34 +01:00
2023-02-06 22:33:26 +01:00
# Process S-vertices waiting to be scanned.
# This loop runs through O(n) iterations per stage.
for x in self.scan_queue:
2023-02-04 21:06:34 +01:00
2023-02-07 11:23:42 +01:00
# Double-check that "x" is an S-vertex.
bx = self.vertex_set_node[x].find()
assert bx.label == _LABEL_S
2023-02-04 21:06:34 +01:00
2023-02-07 11:23:42 +01:00
# Scan the edges that are incident on "x".
# This loop runs through O(m) iterations per stage.
2023-02-07 11:23:42 +01:00
for e in adjacent_edges[x]:
(p, q, _w) = edges[e]
2023-02-07 16:15:59 +01:00
y = p if p != x else q
2023-02-04 21:06:34 +01:00
2023-02-07 11:23:42 +01:00
# Consider the edge between vertices "x" and "y".
# Update delta2 or delta3 tracking accordingly.
#
# We don't actually use the edge right now to extend
# the alternating tree or create a blossom or alternating path.
# If appropriate, insert this edge into delta2 or delta3
# tracking.
# Insert this edge into delta2 or delta3 tracking
2023-02-06 22:33:26 +01:00
# Try to pull this edge into an alternating tree.
2023-02-04 21:06:34 +01:00
2023-02-06 22:33:26 +01:00
# Ignore edges that are internal to a blossom.
by = self.vertex_set_node[y].find()
if bx is by:
2023-02-06 22:33:26 +01:00
continue
2023-02-04 21:06:34 +01:00
if by.label == _LABEL_S:
2024-06-29 22:06:28 +02:00
self.delta3_add_edge(e)
else:
2024-06-29 22:06:28 +02:00
self.delta2_add_edge(e, y, by)
2023-02-06 22:33:26 +01:00
self.scan_queue.clear()
#
# Delta steps:
#
def calc_dual_delta_step(
self
) -> tuple[int, float, int, Optional[_NonTrivialBlossom]]:
2023-02-06 22:33:26 +01:00
"""Calculate a delta step in the dual LPP problem.
2023-02-04 21:06:34 +01:00
2023-02-06 22:33:26 +01:00
This function returns the minimum of the 4 types of delta values,
and the type of delta which obtain the minimum, and the edge or
blossom that produces the minimum delta, if applicable.
2023-02-04 21:06:34 +01:00
2023-02-06 22:33:26 +01:00
The returned value is 2 times the actual delta value.
Multiplication by 2 ensures that the result is an integer if all edge
weights are integers.
2023-02-04 21:06:34 +01:00
This function takes time O((1 + k) * log(n)),
where "k" is the number of intra-blossom edges removed from
the delta3 queue.
2024-06-30 21:38:12 +02:00
At most O(n) delta steps can occur during a stage.
Each edge can be inserted into the delta3 queue at most once per stage.
Therefore, this function takes total time O((n + m) * log(n))
per stage.
2023-02-04 21:06:34 +01:00
2023-02-06 22:33:26 +01:00
Returns:
Tuple (delta_type, delta_2x, delta_edge, delta_blossom).
"""
delta_edge = -1
delta_blossom: Optional[_NonTrivialBlossom] = None
2023-02-06 22:33:26 +01:00
# Compute delta1: minimum dual variable of any S-vertex.
# All unmatched vertices have the same dual value, and this is
# the minimum value among all S-vertices.
2023-02-06 22:33:26 +01:00
delta_type = 1
delta_2x = self.start_vertex_dual_2x - self.delta_sum_2x
2023-02-06 22:33:26 +01:00
# Compute delta2: minimum slack of any edge between an S-vertex and
# an unlabeled vertex.
# This takes time O(log(n)).
2024-06-29 22:06:28 +02:00
(e, slack) = self.delta2_get_min_edge()
if (e != -1) and (slack <= delta_2x):
delta_type = 2
delta_2x = slack
delta_edge = e
2023-02-06 22:33:26 +01:00
# Compute delta3: half minimum slack of any edge between two top-level
# S-blossoms.
2024-06-29 22:06:28 +02:00
# This takes total time O(m * log(n)) per stage.
(e, slack) = self.delta3_get_min_edge()
if (e != -1) and (slack <= delta_2x):
delta_type = 3
delta_2x = slack
delta_edge = e
2023-02-06 22:33:26 +01:00
# Compute delta4: half minimum dual variable of a top-level T-blossom.
# This takes time O(log(n)).
if not self.delta4_queue.empty():
blossom = self.delta4_queue.find_min().data
assert blossom.label == _LABEL_T
assert blossom.parent is None
2024-05-26 11:32:02 +02:00
blossom_dual = blossom.dual_var - self.delta_sum_2x
if blossom_dual <= delta_2x:
delta_type = 4
2024-05-26 11:32:02 +02:00
delta_2x = blossom_dual
delta_blossom = blossom
2023-02-04 21:06:34 +01:00
2023-02-06 22:33:26 +01:00
return (delta_type, delta_2x, delta_edge, delta_blossom)
2023-02-04 21:06:34 +01:00
#
# Main algorithm:
#
2023-02-06 22:33:26 +01:00
def start(self) -> None:
"""Mark each vertex as the node of an alternating tree.
Assign label S to all vertices and add them to the scan queue.
2024-06-30 21:38:12 +02:00
This function takes time O(n + m).
It is called once, at the beginning of the algorithm.
"""
for x in range(self.graph.num_vertex):
assert self.vertex_mate[x] == -1
bx = self.vertex_set_node[x].find()
assert bx.base_vertex == x
# Assign label S.
self.assign_blossom_label_s(bx)
# Mark blossom as the root of an alternating tree.
bx.tree_edge = None
bx.tree_blossoms = {bx}
2023-02-06 22:33:26 +01:00
def run_stage(self) -> bool:
"""Run one stage of the matching algorithm.
The stage searches a maximum-weight augmenting path.
If this path is found, it is used to augment the matching,
thereby increasing the number of matched edges by 1.
If no such path is found, the matching must already be optimal.
This function takes time O((n + m) * log(n)).
2023-02-06 22:33:26 +01:00
Returns:
True if the matching was successfully augmented.
False if no further improvement is possible.
"""
2023-02-05 19:33:13 +01:00
2023-02-06 22:33:26 +01:00
# Each pass through the following loop is a "substage".
# The substage tries to find an augmenting path.
# If an augmenting path is found, we augment the matching and end
# the stage. Otherwise we update the dual LPP problem and enter the
# next substage, or stop if no further improvement is possible.
#
# This loop runs through at most O(n) iterations per stage.
while True:
2023-02-04 21:06:34 +01:00
2024-07-06 16:13:34 +02:00
# self._check_alternating_tree_consistency() # TODO -- remove this
# Consider the incident edges of newly labeled S-vertices.
self.scan_new_s_vertices()
2023-02-04 21:06:34 +01:00
2023-02-06 22:33:26 +01:00
# Calculate delta step in the dual LPP problem.
(delta_type, delta_2x, delta_edge, delta_blossom
) = self.calc_dual_delta_step()
2023-02-06 22:33:26 +01:00
2024-05-26 20:55:54 +02:00
# Update the running sum of delta steps.
# This implicitly updates the dual variables as needed, because
# the running delta sum is taken into account when calculating
# dual values.
self.delta_sum_2x += delta_2x
2023-02-06 22:33:26 +01:00
if delta_type == 2:
# Use the edge from S-vertex to unlabeled vertex that got
# unlocked through the delta update.
2023-02-07 11:23:42 +01:00
(x, y, _w) = self.graph.edges[delta_edge]
if self.vertex_set_node[x].find().label != _LABEL_S:
2023-02-07 11:23:42 +01:00
(x, y) = (y, x)
self.extend_tree_s_to_t(x, y)
2023-02-06 22:33:26 +01:00
elif delta_type == 3:
# Use the S-to-S edge that got unlocked by the delta update.
# This reveals either a new blossom or an augmenting path.
2023-02-07 11:23:42 +01:00
(x, y, _w) = self.graph.edges[delta_edge]
if self.add_s_to_s_edge(x, y):
# Matching was augmented. End the stage.
return True
2023-02-06 22:33:26 +01:00
elif delta_type == 4:
# Expand the T-blossom that reached dual value 0 through
# the delta update.
assert delta_blossom is not None
2023-02-06 22:33:26 +01:00
self.expand_t_blossom(delta_blossom)
2023-02-04 21:06:34 +01:00
2023-02-06 22:33:26 +01:00
else:
# No further improvement possible. End the algorithm.
2023-02-06 22:33:26 +01:00
assert delta_type == 1
return False
2023-02-06 22:33:26 +01:00
def cleanup(self) -> None:
"""Remove all alternating trees and mark all blossoms as unlabeled.
2023-02-04 21:06:34 +01:00
Also applies delayed updates to dual variables.
Also resets tracking of least-slack edges.
2024-06-30 21:38:12 +02:00
This function takes time O((n + m) * log(n)).
2024-06-25 21:01:52 +02:00
It is called once, at the end of the algorithm.
"""
2024-06-25 21:01:52 +02:00
assert not self.scan_queue
for blossom in itertools.chain(self.trivial_blossom,
self.nontrivial_blossom):
# Remove blossom label.
if (blossom.parent is None) and (blossom.label != _LABEL_NONE):
self.reset_blossom_label(blossom)
assert blossom.label == _LABEL_NONE
# Remove blossom from alternating tree.
blossom.tree_edge = None
blossom.tree_blossoms = None
# Unwind lazy delta updates to vertex dual variables.
if blossom.vertex_dual_offset != 0:
for x in blossom.vertices():
self.vertex_dual_2x[x] += blossom.vertex_dual_offset
blossom.vertex_dual_offset = 0
assert self.delta2_queue.empty()
assert self.delta3_queue.empty()
assert self.delta4_queue.empty()
2023-02-04 21:06:34 +01:00
def _verify_blossom_edges(
ctx: _MatchingContext,
blossom: _NonTrivialBlossom,
edge_slack_2x: list[float]
) -> None:
"""Descend down the blossom tree to find edges that are contained
in blossoms.
Adjust the slack of all contained edges to account for the dual variables
of its containing blossoms.
On the way down, keep track of the sum of dual variables of
the containing blossoms.
On the way up, keep track of the total number of matched edges
in the subblossoms. Then check that all blossoms with non-zero
dual variable are "full".
Raises:
2023-02-12 21:18:57 +01:00
MatchingError: If a blossom with non-zero dual is not full.
"""
num_vertex = ctx.graph.num_vertex
# For each vertex "x",
# "vertex_depth[x]" is the depth of the smallest blossom on
# the current descent path that contains "x".
vertex_depth: list[int] = num_vertex * [0]
# Keep track of the sum of blossom duals at each depth along
# the current descent path.
path_sum_dual: list[float] = [0]
# Keep track of the number of matched edges at each depth along
# the current descent path.
path_num_matched: list[int] = [0]
# Use an explicit stack to avoid deep recursion.
stack: list[tuple[_NonTrivialBlossom, int]] = [(blossom, -1)]
while stack:
(blossom, p) = stack[-1]
depth = len(stack)
if p == -1:
# We just entered this sub-blossom.
# Update the depth of all vertices in this sub-blossom.
for x in blossom.vertices():
vertex_depth[x] = depth
# Calculate the sub of blossoms at the current depth.
path_sum_dual.append(path_sum_dual[-1] + blossom.dual_var)
# Initialize the number of matched edges at the current depth.
path_num_matched.append(0)
p += 1
if p < len(blossom.subblossoms):
# Update the sub-blossom pointer at the current level.
stack[-1] = (blossom, p + 1)
# Examine the next sub-blossom at the current level.
sub = blossom.subblossoms[p]
if isinstance(sub, _NonTrivialBlossom):
# Prepare to descent into the selected sub-blossom and
# scan it recursively.
stack.append((sub, -1))
else:
# Handle this trivial sub-blossom.
# Scan its adjacent edges and find the smallest blossom
# that contains each edge.
for e in ctx.graph.adjacent_edges[sub.base_vertex]:
(x, y, _w) = ctx.graph.edges[e]
# Only process edges that are ordered out from this
# sub-blossom. This ensures that we process each edge in
# the blossom only once.
if x == sub.base_vertex:
edge_depth = vertex_depth[y]
if edge_depth > 0:
# This edge is contained in an ancestor blossom.
# Update its slack.
edge_slack_2x[e] += 2 * path_sum_dual[edge_depth]
# Update the number of matched edges in ancestor.
if ctx.vertex_mate[x] == y:
path_num_matched[edge_depth] += 1
else:
# We are now leaving the current sub-blossom.
# Count the number of vertices inside this blossom.
blossom_vertices = blossom.vertices()
blossom_num_vertex = len(blossom_vertices)
# Check that all blossoms are "full".
# A blossom is full if all except one of its vertices are
# matched to another vertex in the blossom.
blossom_num_matched = path_num_matched[depth]
if blossom_num_vertex != 2 * blossom_num_matched + 1:
raise MatchingError(
"Verification failed: blossom non-full"
f" dual={blossom.dual_var}"
f" nvertex={blossom_num_vertex}"
f" nmatched={blossom_num_matched}")
# Update the number of matched edges in the parent blossom to
# take into account the matched edges in this blossom.
path_num_matched[depth - 1] += path_num_matched[depth]
# Revert the depth of the vertices in this sub-blossom.
for x in blossom_vertices:
vertex_depth[x] = depth - 1
# Trim the descending path.
path_sum_dual.pop()
path_num_matched.pop()
# Remove the current blossom from the stack.
# We thus continue our scan of the parent blossom.
stack.pop()
2023-02-06 22:33:26 +01:00
def _verify_optimum(ctx: _MatchingContext) -> None:
2023-02-04 21:06:34 +01:00
"""Verify that the optimum solution has been found.
This function takes time O(n**2).
2023-02-04 21:06:34 +01:00
Raises:
2023-02-12 21:18:57 +01:00
MatchingError: If the solution is not optimal.
2023-02-04 21:06:34 +01:00
"""
2023-02-06 22:33:26 +01:00
num_vertex = ctx.graph.num_vertex
num_edge = len(ctx.graph.edges)
2023-02-12 21:18:57 +01:00
# Check that each matched edge actually exists in the graph.
num_matched_vertex = 0
2023-02-07 11:23:42 +01:00
for x in range(num_vertex):
2023-02-12 21:18:57 +01:00
y = ctx.vertex_mate[x]
if y != -1:
if ctx.vertex_mate[y] != x:
raise MatchingError(
"Verification failed:"
f" asymmetric match of vertex {x} and {y}")
num_matched_vertex += 1
num_matched_edge = 0
2023-02-07 11:23:42 +01:00
for (x, y, _w) in ctx.graph.edges:
if ctx.vertex_mate[x] == y:
num_matched_edge += 1
2023-02-12 21:18:57 +01:00
if num_matched_vertex != 2 * num_matched_edge:
raise MatchingError(
f"Verification failed: {num_matched_vertex} matched vertices"
f" inconsistent with {num_matched_edge} matched edges")
2023-02-04 21:06:34 +01:00
# Check that all dual variables are non-negative.
for x in range(num_vertex):
if ctx.vertex_dual_2x[x] < 0:
raise MatchingError(
"Verification failed:"
f" vertex {x} has negative dual {ctx.vertex_dual_2x[x]/2}")
for blossom in ctx.nontrivial_blossom:
2023-02-12 21:18:57 +01:00
if blossom.dual_var < 0:
raise MatchingError("Verification failed:"
f" negative blossom dual {blossom.dual_var}")
2023-02-04 21:06:34 +01:00
# Check that all unmatched vertices have zero dual.
for x in range(num_vertex):
if ctx.vertex_mate[x] == -1 and ctx.vertex_dual_2x[x] != 0:
raise MatchingError(
f"Verification failed: Unmatched vertex {x}"
f" has non-zero dual {ctx.vertex_dual_2x[x]/2}")
# Calculate the slack of each edge.
# A correction will be needed for edges inside blossoms.
edge_slack_2x: list[float] = [
ctx.vertex_dual_2x[x] + ctx.vertex_dual_2x[y] - 2 * w
for (x, y, w) in ctx.graph.edges]
2023-02-04 21:06:34 +01:00
# Descend down each top-level blossom.
# Adjust edge slacks to account for the duals of its containing blossoms.
# And check that all blossoms are full.
# This takes total time O(n**2).
for blossom in ctx.nontrivial_blossom:
if blossom.parent is None:
_verify_blossom_edges(ctx, blossom, edge_slack_2x)
# We now know the correct slack of each edge.
# Check that all edges have non-negative slack.
2023-02-12 21:18:57 +01:00
min_edge_slack = min(edge_slack_2x)
if min_edge_slack < 0:
raise MatchingError(
f"Verification failed: negative edge slack {min_edge_slack/2}")
# Check that all matched edges have zero slack.
for e in range(num_edge):
(x, y, _w) = ctx.graph.edges[e]
2023-02-12 21:18:57 +01:00
if ctx.vertex_mate[x] == y and edge_slack_2x[e] != 0:
raise MatchingError(
"Verification failed:"
f" matched edge ({x}, {y}) has slack {edge_slack_2x[e]/2}")
2023-02-04 21:06:34 +01:00
# Optimum solution confirmed.