deps: replace hopcroftkarp with scipy.sparse.csgraph.maximum_bipartite_matching

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:

pyproject.toml

@@ -24,11 +24,11 @@ maintainers = [
 
 dependencies = [
   "deprecated",
-  "hopcroftkarp",
   "joblib",
   "matplotlib",
   "numpy",
   "scikit-learn",
+  "scipy",
 ]
 
 classifiers = [

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):