From eceb5842f734641acc39fe2ddaa070180b543087 Mon Sep 17 00:00:00 2001 From: Sushovan Majhi Date: Sun, 9 Aug 2026 15:16:12 -0400 Subject: [PATCH] deps: replace hopcroftkarp with scipy.sparse.csgraph Closes the two changes requested in #106. `hopcroftkarp` is GPLv3 and was last released 2019-10-11, so no install carrying persim can offer a permissive dependency closure -- and since ripser depends on persim, that propagates across much of the scikit-tda stack. `scipy.sparse.csgraph.maximum_bipartite_matching` is Hopcroft-Karp, and scipy is already present in every working install via scikit-learn, though it was undeclared. So this removes a dependency and adds none: the scipy entry makes an existing requirement visible. The call site changes shape rather than name. hopcroftkarp took a dict of string-keyed adjacency sets and returned a dict holding both directions of every matched pair, so a perfect matching was tested as len(res) == 2 * n. scipy takes a CSR matrix and returns an array where res[i] is the column matched to row i, or -1 -- so the test becomes np.all(res >= 0), which is equivalent here because D is square. Building the CSR matrix from a boolean array drops non-edges rather than storing them as explicit zeros, which is what maximum_bipartite_matching reads. Bottleneck distances are unchanged: verified identical on the existing suite and on 40 randomly generated diagram pairs. The returned matching may differ where more than one is optimal. The distance is unique; the matching attaining it need not be, and the two routines break ties differently. On 17 of those 40 pairs the matching differs while the distance does not, and the row count can differ too, since a diagonal-to-diagonal row is dropped on the way out. Both remain valid: every point of both diagrams is accounted for exactly once and no pair exceeds the distance. test_matching_is_ optimal_when_several_are pins that contract on a pair where the choice is real, and deliberately does not pin the pairs themselves. --- persim/bottleneck.py | 23 ++++++++++++++--------- pyproject.toml | 2 +- test/test_distances.py | 33 +++++++++++++++++++++++++++++++++ 3 files changed, 48 insertions(+), 10 deletions(-) diff --git a/persim/bottleneck.py b/persim/bottleneck.py index 286e668b..fdf6ebe8 100644 --- a/persim/bottleneck.py +++ b/persim/bottleneck.py @@ -1,7 +1,8 @@ """ Implementation of the bottleneck distance using binary -search and the Hopcroft-Karp algorithm +search and the Hopcroft-Karp algorithm, via +scipy.sparse.csgraph.maximum_bipartite_matching Author: Chris Tralie @@ -10,7 +11,8 @@ import numpy as np from bisect import bisect_left -from hopcroftkarp import HopcroftKarp +from scipy.sparse import csr_matrix +from scipy.sparse.csgraph import maximum_bipartite_matching import warnings __all__ = ["bottleneck"] @@ -100,17 +102,20 @@ def bottleneck(dgm1, dgm2, matching=False): # bottleneck distance ds = np.sort(np.unique(D.flatten())) # [0:-1] # Everything but np.inf bdist = ds[-1] - matching = {} + matching = np.array([], dtype=int) while len(ds) >= 1: idx = 0 if len(ds) > 1: idx = bisect_left(range(ds.size), int(ds.size / 2)) d = ds[idx] - graph = {} - for i in range(D.shape[0]): - graph["{}".format(i)] = {j for j in range(D.shape[1]) if D[i, j] <= d} - res = HopcroftKarp(graph).maximum_matching() - if len(res) == 2 * D.shape[0] and d <= bdist: + # Edges are the pairs within distance d. `maximum_bipartite_matching` + # reads the sparsity structure, and building the CSR matrix from a + # boolean array drops the non-edges rather than storing them as + # explicit zeros. + res = maximum_bipartite_matching(csr_matrix(D <= d), perm_type="column") + # res[i] is the column matched to row i, or -1 if row i is unmatched. + # D is square, so every row being matched is a perfect matching. + if np.all(res >= 0) and d <= bdist: bdist = d matching = res ds = ds[0:idx] @@ -120,7 +125,7 @@ def bottleneck(dgm1, dgm2, matching=False): if return_matching: matchidx = [] for i in range(M + N): - j = matching["{}".format(i)] + j = matching[i] d = D[i, j] if i < M: if j >= N: diff --git a/pyproject.toml b/pyproject.toml index c1e27dc8..abea31af 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -24,11 +24,11 @@ maintainers = [ dependencies = [ "deprecated", - "hopcroftkarp", "joblib", "matplotlib", "numpy", "scikit-learn", + "scipy", ] classifiers = [ diff --git a/test/test_distances.py b/test/test_distances.py index 94346917..168e501d 100644 --- a/test/test_distances.py +++ b/test/test_distances.py @@ -110,6 +110,39 @@ def test_one_diagonal(self): assert dist == 5.0 assert returned_matching.shape[1] == 3 + def test_matching_is_optimal_when_several_are(self): + # The bottleneck distance is unique; the optimal matching attaining it + # is not. Where more than one attains it, which one is returned is a + # property of the matching routine rather than of the diagrams, so this + # pins the contract every optimal matching satisfies and deliberately + # does not pin the pairs themselves. + # + # On this pair the choice is real rather than hypothetical: the cost is + # set by dgm1's third point going to the diagonal, which leaves dgm1's + # first point free to pair with dgm2's only point or to go to the + # diagonal as well. Both are optimal, and they differ in how many rows + # the matching has -- a diagonal-to-diagonal row is dropped on the way + # out, so asserting a row count here would assert an implementation. + dgm1 = np.array([[0.9726, 1.4526], [0.8899, 1.1223], [0.8224, 1.6243]]) + dgm2 = np.array([[0.9235, 1.1897]]) + + dist, matching = bottleneck(dgm1, dgm2, matching=True) + + # Set by dgm1[2] -> diagonal, at 0.5 * (1.6243 - 0.8224). + assert dist == pytest.approx(0.40095, abs=1e-9) + + # Every point of both diagrams is accounted for exactly once, whether + # it was paired across or sent to the diagonal. + left = np.sort(matching[matching[:, 0] >= 0, 0]) + right = np.sort(matching[matching[:, 1] >= 0, 1]) + assert left.tolist() == list(range(dgm1.shape[0])) + assert right.tolist() == list(range(dgm2.shape[0])) + + # No pair in an optimal matching may cost more than the distance the + # matching attains, and at least one must attain it. + assert np.all(matching[:, 2] <= dist + 1e-12) + assert np.max(matching[:, 2]) == pytest.approx(dist) + class TestWasserstein: def test_single(self):