1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
|
# This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT.
# See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details.
# Author(s): Mathieu Carrière
#
# Copyright (C) 2018-2019 Inria
#
# Modification(s):
# - YYYY/MM Author: Description of the modification
import numpy as np
from sklearn.base import BaseEstimator, TransformerMixin
from sklearn.metrics import pairwise_distances
try:
from .. import bottleneck_distance
USE_GUDHI = True
except ImportError:
USE_GUDHI = False
print("Gudhi built without CGAL: BottleneckDistance will return a null matrix")
#############################################
# Metrics ###################################
#############################################
class SlicedWassersteinDistance(BaseEstimator, TransformerMixin):
"""
This is a class for computing the sliced Wasserstein distance matrix from a list of persistence diagrams. The Sliced Wasserstein distance is computed by projecting the persistence diagrams onto lines, comparing the projections with the 1-norm, and finally integrating over all possible lines. See http://proceedings.mlr.press/v70/carriere17a.html for more details.
"""
def __init__(self, num_directions=10):
"""
Constructor for the SlicedWassersteinDistance class.
Parameters:
num_directions (int): number of lines evenly sampled from [-pi/2,pi/2] in order to approximate and speed up the distance computation (default 10).
"""
self.num_directions = num_directions
thetas = np.linspace(-np.pi/2, np.pi/2, num=self.num_directions+1)[np.newaxis,:-1]
self.lines_ = np.concatenate([np.cos(thetas), np.sin(thetas)], axis=0)
def fit(self, X, y=None):
"""
Fit the SlicedWassersteinDistance class on a list of persistence diagrams: persistence diagrams are projected onto the different lines. The diagrams themselves and their projections are then stored in numpy arrays, called **diagrams_** and **approx_diag_**.
Parameters:
X (list of n x 2 numpy arrays): input persistence diagrams.
y (n x 1 array): persistence diagram labels (unused).
"""
self.diagrams_ = X
self.approx_ = [np.matmul(X[i], self.lines_) for i in range(len(X))]
diag_proj = (1./2) * np.ones((2,2))
self.approx_diag_ = [np.matmul(np.matmul(X[i], diag_proj), self.lines_) for i in range(len(X))]
return self
def transform(self, X):
"""
Compute all sliced Wasserstein distances between the persistence diagrams that were stored after calling the fit() method, and a given list of (possibly different) persistence diagrams.
Parameters:
X (list of n x 2 numpy arrays): input persistence diagrams.
Returns:
numpy array of shape (number of diagrams in **diagrams**) x (number of diagrams in X): matrix of pairwise sliced Wasserstein distances.
"""
Xfit = np.zeros((len(X), len(self.approx_)))
if len(self.diagrams_) == len(X) and np.all([np.array_equal(self.diagrams_[i], X[i]) for i in range(len(X))]):
for i in range(len(self.approx_)):
for j in range(i+1, len(self.approx_)):
A = np.sort(np.concatenate([self.approx_[i], self.approx_diag_[j]], axis=0), axis=0)
B = np.sort(np.concatenate([self.approx_[j], self.approx_diag_[i]], axis=0), axis=0)
L1 = np.sum(np.abs(A-B), axis=0)
Xfit[i,j] = np.mean(L1)
Xfit[j,i] = Xfit[i,j]
else:
diag_proj = (1./2) * np.ones((2,2))
approx = [np.matmul(X[i], self.lines_) for i in range(len(X))]
approx_diag = [np.matmul(np.matmul(X[i], diag_proj), self.lines_) for i in range(len(X))]
for i in range(len(approx)):
for j in range(len(self.approx_)):
A = np.sort(np.concatenate([approx[i], self.approx_diag_[j]], axis=0), axis=0)
B = np.sort(np.concatenate([self.approx_[j], approx_diag[i]], axis=0), axis=0)
L1 = np.sum(np.abs(A-B), axis=0)
Xfit[i,j] = np.mean(L1)
return Xfit
class BottleneckDistance(BaseEstimator, TransformerMixin):
"""
This is a class for computing the bottleneck distance matrix from a list of persistence diagrams.
"""
def __init__(self, epsilon=None):
"""
Constructor for the BottleneckDistance class.
Parameters:
epsilon (double): absolute (additive) error tolerated on the distance (default is the smallest positive float), see :func:`gudhi.bottleneck_distance`.
"""
self.epsilon = epsilon
def fit(self, X, y=None):
"""
Fit the BottleneckDistance class on a list of persistence diagrams: persistence diagrams are stored in a numpy array called **diagrams**.
Parameters:
X (list of n x 2 numpy arrays): input persistence diagrams.
y (n x 1 array): persistence diagram labels (unused).
"""
self.diagrams_ = X
return self
def transform(self, X):
"""
Compute all bottleneck distances between the persistence diagrams that were stored after calling the fit() method, and a given list of (possibly different) persistence diagrams.
Parameters:
X (list of n x 2 numpy arrays): input persistence diagrams.
Returns:
numpy array of shape (number of diagrams in **diagrams**) x (number of diagrams in X): matrix of pairwise bottleneck distances.
"""
num_diag1 = len(X)
#if len(self.diagrams_) == len(X) and np.all([np.array_equal(self.diagrams_[i], X[i]) for i in range(len(X))]):
if X is self.diagrams_:
matrix = np.zeros((num_diag1, num_diag1))
if USE_GUDHI:
for i in range(num_diag1):
for j in range(i+1, num_diag1):
matrix[i,j] = bottleneck_distance(X[i], X[j], self.epsilon)
matrix[j,i] = matrix[i,j]
else:
print("Gudhi built without CGAL: returning a null matrix")
else:
num_diag2 = len(self.diagrams_)
matrix = np.zeros((num_diag1, num_diag2))
if USE_GUDHI:
for i in range(num_diag1):
for j in range(num_diag2):
matrix[i,j] = bottleneck_distance(X[i], self.diagrams_[j], self.epsilon)
else:
print("Gudhi built without CGAL: returning a null matrix")
Xfit = matrix
return Xfit
class PersistenceFisherDistance(BaseEstimator, TransformerMixin):
"""
This is a class for computing the persistence Fisher distance matrix from a list of persistence diagrams. The persistence Fisher distance is obtained by computing the original Fisher distance between the probability distributions associated to the persistence diagrams given by convolving them with a Gaussian kernel. See http://papers.nips.cc/paper/8205-persistence-fisher-kernel-a-riemannian-manifold-kernel-for-persistence-diagrams for more details.
"""
def __init__(self, bandwidth=1., kernel_approx=None):
"""
Constructor for the PersistenceFisherDistance class.
Parameters:
bandwidth (double): bandwidth of the Gaussian kernel used to turn persistence diagrams into probability distributions (default 1.).
kernel_approx (class): kernel approximation class used to speed up computation (default None). Common kernel approximations classes can be found in the scikit-learn library (such as RBFSampler for instance).
"""
self.bandwidth, self.kernel_approx = bandwidth, kernel_approx
def fit(self, X, y=None):
"""
Fit the PersistenceFisherDistance class on a list of persistence diagrams: persistence diagrams are stored in a numpy array called **diagrams** and the kernel approximation class (if not None) is applied on them.
Parameters:
X (list of n x 2 numpy arrays): input persistence diagrams.
y (n x 1 array): persistence diagram labels (unused).
"""
self.diagrams_ = X
projection = (1./2) * np.ones((2,2))
self.diagonal_projections_ = [np.matmul(X[i], projection) for i in range(len(X))]
if self.kernel_approx is not None:
self.approx_ = [self.kernel_approx.transform(X[i]) for i in range(len(X))]
self.approx_diagonal_ = [self.kernel_approx.transform(self.diagonal_projections_[i]) for i in range(len(X))]
return self
def transform(self, X):
"""
Compute all persistence Fisher distances between the persistence diagrams that were stored after calling the fit() method, and a given list of (possibly different) persistence diagrams.
Parameters:
X (list of n x 2 numpy arrays): input persistence diagrams.
Returns:
numpy array of shape (number of diagrams in **diagrams**) x (number of diagrams in X): matrix of pairwise persistence Fisher distances.
"""
Xfit = np.zeros((len(X), len(self.diagrams_)))
if len(self.diagrams_) == len(X) and np.all([np.array_equal(self.diagrams_[i], X[i]) for i in range(len(X))]):
for i in range(len(self.diagrams_)):
for j in range(i+1, len(self.diagrams_)):
if self.kernel_approx is not None:
Z = np.concatenate([self.approx_[i], self.approx_diagonal_[i], self.approx_[j], self.approx_diagonal_[j]], axis=0)
U, V = np.sum(np.concatenate([self.approx_[i], self.approx_diagonal_[j]], axis=0), axis=0), np.sum(np.concatenate([self.approx_[j], self.approx_diagonal_[i]], axis=0), axis=0)
vectori, vectorj = np.abs(np.matmul(Z, U.T)), np.abs(np.matmul(Z, V.T))
vectori_sum, vectorj_sum = np.sum(vectori), np.sum(vectorj)
if vectori_sum != 0:
vectori = vectori/vectori_sum
if vectorj_sum != 0:
vectorj = vectorj/vectorj_sum
Xfit[i,j] = np.arccos( min(np.dot(np.sqrt(vectori), np.sqrt(vectorj)), 1.) )
Xfit[j,i] = Xfit[i,j]
else:
Z = np.concatenate([self.diagrams_[i], self.diagonal_projections_[i], self.diagrams_[j], self.diagonal_projections_[j]], axis=0)
U, V = np.concatenate([self.diagrams_[i], self.diagonal_projections_[j]], axis=0), np.concatenate([self.diagrams_[j], self.diagonal_projections_[i]], axis=0)
vectori = np.sum(np.exp(-np.square(pairwise_distances(Z,U))/(2 * np.square(self.bandwidth)))/(self.bandwidth * np.sqrt(2*np.pi)), axis=1)
vectorj = np.sum(np.exp(-np.square(pairwise_distances(Z,V))/(2 * np.square(self.bandwidth)))/(self.bandwidth * np.sqrt(2*np.pi)), axis=1)
vectori_sum, vectorj_sum = np.sum(vectori), np.sum(vectorj)
if vectori_sum != 0:
vectori = vectori/vectori_sum
if vectorj_sum != 0:
vectorj = vectorj/vectorj_sum
Xfit[i,j] = np.arccos( min(np.dot(np.sqrt(vectori), np.sqrt(vectorj)), 1.) )
Xfit[j,i] = Xfit[i,j]
else:
projection = (1./2) * np.ones((2,2))
diagonal_projections = [np.matmul(X[i], projection) for i in range(len(X))]
if self.kernel_approx is not None:
approx = [self.kernel_approx.transform(X[i]) for i in range(len(X))]
approx_diagonal = [self.kernel_approx.transform(diagonal_projections[i]) for i in range(len(X))]
for i in range(len(X)):
for j in range(len(self.diagrams_)):
if self.kernel_approx is not None:
Z = np.concatenate([approx[i], approx_diagonal[i], self.approx_[j], self.approx_diagonal_[j]], axis=0)
U, V = np.sum(np.concatenate([approx[i], self.approx_diagonal_[j]], axis=0), axis=0), np.sum(np.concatenate([self.approx_[j], approx_diagonal[i]], axis=0), axis=0)
vectori, vectorj = np.abs(np.matmul(Z, U.T)), np.abs(np.matmul(Z, V.T))
vectori_sum, vectorj_sum = np.sum(vectori), np.sum(vectorj)
if vectori_sum != 0:
vectori = vectori/vectori_sum
if vectorj_sum != 0:
vectorj = vectorj/vectorj_sum
Xfit[i,j] = np.arccos( min(np.dot(np.sqrt(vectori), np.sqrt(vectorj)), 1.) )
else:
Z = np.concatenate([X[i], diagonal_projections[i], self.diagrams_[j], self.diagonal_projections_[j]], axis=0)
U, V = np.concatenate([X[i], self.diagonal_projections_[j]], axis=0), np.concatenate([self.diagrams_[j], diagonal_projections[i]], axis=0)
vectori = np.sum(np.exp(-np.square(pairwise_distances(Z,U))/(2 * np.square(self.bandwidth)))/(self.bandwidth * np.sqrt(2*np.pi)), axis=1)
vectorj = np.sum(np.exp(-np.square(pairwise_distances(Z,V))/(2 * np.square(self.bandwidth)))/(self.bandwidth * np.sqrt(2*np.pi)), axis=1)
vectori_sum, vectorj_sum = np.sum(vectori), np.sum(vectorj)
if vectori_sum != 0:
vectori = vectori/vectori_sum
if vectorj_sum != 0:
vectorj = vectorj/vectorj_sum
Xfit[i,j] = np.arccos( min(np.dot(np.sqrt(vectori), np.sqrt(vectorj)), 1.) )
return Xfit
|
