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