scikit-tda / persim
 1 ```""" ``` 2 3 ``` Implementation of the bottleneck distance using binary ``` 4 ``` search and the Hopcroft-Karp algorithm ``` 5 6 ``` Author: Chris Tralie ``` 7 8 ```""" ``` 9 10 3 ```import numpy as np ``` 11 12 3 ```from bisect import bisect_left ``` 13 3 ```from hopcroftkarp import HopcroftKarp ``` 14 3 ```import warnings ``` 15 16 3 ```__all__ = ["bottleneck"] ``` 17 18 19 3 ```def bottleneck(dgm1, dgm2, matching=False): ``` 20 ``` """ ``` 21 ``` Perform the Bottleneck distance matching between persistence diagrams. ``` 22 ``` Assumes first two columns of S and T are the coordinates of the persistence ``` 23 ``` points, but allows for other coordinate columns (which are ignored in ``` 24 ``` diagonal matching). ``` 25 26 ``` See the `distances` notebook for an example of how to use this. ``` 27 28 ``` Parameters ``` 29 ``` ----------- ``` 30 ``` dgm1: Mx(>=2) ``` 31 ``` array of birth/death pairs for PD 1 ``` 32 ``` dgm2: Nx(>=2) ``` 33 ``` array of birth/death paris for PD 2 ``` 34 ``` matching: bool, default False ``` 35 ``` if True, return matching infromation and cross-similarity matrix ``` 36 37 ``` Returns ``` 38 ``` -------- ``` 39 40 ``` d: float ``` 41 ``` bottleneck distance between dgm1 and dgm2 ``` 42 ``` (matching, D): Only returns if `matching=True` ``` 43 ``` (tuples of matched indices, (N+M)x(N+M) cross-similarity matrix) ``` 44 ``` """ ``` 45 46 3 ``` return_matching = matching ``` 47 48 3 ``` S = np.array(dgm1) ``` 49 3 ``` M = min(S.shape[0], S.size) ``` 50 3 ``` if S.size > 0: ``` 51 3 ``` S = S[np.isfinite(S[:, 1]), :] ``` 52 3 ``` if S.shape[0] < M: ``` 53 3 ``` warnings.warn( ``` 54 ``` "dgm1 has points with non-finite death times;"+ ``` 55 ``` "ignoring those points" ``` 56 ``` ) ``` 57 3 ``` M = S.shape[0] ``` 58 3 ``` T = np.array(dgm2) ``` 59 3 ``` N = min(T.shape[0], T.size) ``` 60 3 ``` if T.size > 0: ``` 61 3 ``` T = T[np.isfinite(T[:, 1]), :] ``` 62 3 ``` if T.shape[0] < N: ``` 63 3 ``` warnings.warn( ``` 64 ``` "dgm2 has points with non-finite death times;"+ ``` 65 ``` "ignoring those points" ``` 66 ``` ) ``` 67 3 ``` N = T.shape[0] ``` 68 69 3 ``` if M == 0: ``` 70 3 ``` S = np.array([[0, 0]]) ``` 71 3 ``` M = 1 ``` 72 3 ``` if N == 0: ``` 73 3 ``` T = np.array([[0, 0]]) ``` 74 3 ``` N = 1 ``` 75 76 ``` # Step 1: Compute CSM between S and T, including points on diagonal ``` 77 ``` # L Infinity distance ``` 78 3 ``` Sb, Sd = S[:, 0], S[:, 1] ``` 79 3 ``` Tb, Td = T[:, 0], T[:, 1] ``` 80 3 ``` D1 = np.abs(Sb[:, None] - Tb[None, :]) ``` 81 3 ``` D2 = np.abs(Sd[:, None] - Td[None, :]) ``` 82 3 ``` DUL = np.maximum(D1, D2) ``` 83 84 ``` # Put diagonal elements into the matrix, being mindful that Linfinity ``` 85 ``` # balls meet the diagonal line at a diamond vertex ``` 86 3 ``` D = np.zeros((M + N, M + N)) ``` 87 3 ``` D[0:M, 0:N] = DUL ``` 88 3 ``` UR = np.max(D) * np.ones((M, M)) ``` 89 3 ``` np.fill_diagonal(UR, 0.5 * (S[:, 1] - S[:, 0])) ``` 90 3 ``` D[0:M, N::] = UR ``` 91 3 ``` UL = np.max(D) * np.ones((N, N)) ``` 92 3 ``` np.fill_diagonal(UL, 0.5 * (T[:, 1] - T[:, 0])) ``` 93 3 ``` D[M::, 0:N] = UL ``` 94 95 ``` # Step 2: Perform a binary search + Hopcroft Karp to find the ``` 96 ``` # bottleneck distance ``` 97 3 ``` M = D.shape[0] ``` 98 3 ``` ds = np.sort(np.unique(D.flatten())) ``` 99 3 ``` bdist = ds[-1] ``` 100 3 ``` matching = {} ``` 101 3 ``` while len(ds) >= 1: ``` 102 3 ``` idx = 0 ``` 103 3 ``` if len(ds) > 1: ``` 104 3 ``` idx = bisect_left(range(ds.size), int(ds.size / 2)) ``` 105 3 ``` d = ds[idx] ``` 106 3 ``` graph = {} ``` 107 3 ``` for i in range(M): ``` 108 3 ``` graph["%s" % i] = {j for j in range(M) if D[i, j] <= d} ``` 109 3 ``` res = HopcroftKarp(graph).maximum_matching() ``` 110 3 ``` if len(res) == 2 * M and d <= bdist: ``` 111 3 ``` bdist = d ``` 112 3 ``` matching = res ``` 113 3 ``` ds = ds[0:idx] ``` 114 ``` else: ``` 115 3 ``` ds = ds[idx + 1::] ``` 116 117 3 ``` if return_matching: ``` 118 3 ``` matchidx = [(i, matching["%i" % i]) for i in range(M)] ``` 119 3 ``` return bdist, (matchidx, D) ``` 120 ``` else: ``` 121 3 ``` return bdist ```

Read our documentation on viewing source code .