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 .

Loading