From a145c7168fdb3f4205cb68870f06fc5cb8e08dea Mon Sep 17 00:00:00 2001
From: MathieuCarriere <mathieu.carriere3@gmail.com>
Date: Fri, 31 Jan 2020 14:49:59 -0500
Subject: factorization of distance and kernel computations

---
 src/python/gudhi/representations/kernel_methods.py | 131 +++++++----
 src/python/gudhi/representations/metrics.py        | 247 +++++++++------------
 2 files changed, 193 insertions(+), 185 deletions(-)

diff --git a/src/python/gudhi/representations/kernel_methods.py b/src/python/gudhi/representations/kernel_methods.py
index bfc83aff..bbbb7c31 100644
--- a/src/python/gudhi/representations/kernel_methods.py
+++ b/src/python/gudhi/representations/kernel_methods.py
@@ -9,13 +9,83 @@
 
 import numpy as np
 from sklearn.base import BaseEstimator, TransformerMixin
-from sklearn.metrics import pairwise_distances
-from .metrics import SlicedWassersteinDistance, PersistenceFisherDistance
+from sklearn.metrics import pairwise_distances, pairwise_kernels
+from .metrics import SlicedWassersteinDistance, PersistenceFisherDistance, sklearn_wrapper, pairwise_persistence_diagram_distances, sliced_wasserstein_distance, persistence_fisher_distance
+from .preprocessing import Padding
 
 #############################################
 # Kernel methods ############################
 #############################################
 
+def persistence_weighted_gaussian_kernel(D1, D2, weight=lambda x: 1, kernel_approx=None, bandwidth=1.):
+    """
+    This is a function for computing the persistence weighted Gaussian kernel value from two persistence diagrams. The persistence weighted Gaussian kernel is computed by convolving the persistence diagram points with weighted Gaussian kernels. See http://proceedings.mlr.press/v48/kusano16.html for more details.
+    :param D1: (n x 2) numpy.array encoding the (finite points of the) first diagram. Must not contain essential points (i.e. with infinite coordinate).
+    :param D2: (m x 2) numpy.array encoding the second diagram.
+    :param bandwidth: bandwidth of the Gaussian kernel with which persistence diagrams will be convolved
+    :param weight: weight function for the persistence diagram points. This function must be defined on 2D points, ie lists or numpy arrays of the form [p_x,p_y].
+    :param kernel_approx: kernel approximation class used to speed up computation. Common kernel approximations classes can be found in the scikit-learn library (such as RBFSampler for instance).
+    :returns: the persistence weighted Gaussian kernel value between persistence diagrams.
+    :rtype: float 
+    """
+    ws1 = np.array([weight(D1[j,:]) for j in range(len(D1))])
+    ws2 = np.array([weight(D2[j,:]) for j in range(len(D2))])
+    if kernel_approx is not None:
+        approx1 = np.sum(np.multiply(ws1[:,np.newaxis], kernel_approx.transform(D1)), axis=0)
+        approx2 = np.sum(np.multiply(ws2[:,np.newaxis], kernel_approx.transform(D2)), axis=0)
+        return (1./(np.sqrt(2*np.pi)*bandwidth)) * np.matmul(approx1, approx2.T)
+    else:
+        W = np.matmul(ws1[:,np.newaxis], ws2[np.newaxis,:])
+        E = (1./(np.sqrt(2*np.pi)*bandwidth)) * np.exp(-np.square(pairwise_distances(D1,D2))/(2*bandwidth*bandwidth))
+        return np.sum(np.multiply(W, E))
+
+def persistence_scale_space_kernel(D1, D2, kernel_approx=None, bandwidth=1.):
+    """
+    This is a function for computing the persistence scale space kernel value from two persistence diagrams. The persistence scale space kernel is computed by adding the symmetric to the diagonal of each point in each persistence diagram, with negative weight, and then convolving the points with a Gaussian kernel. See https://www.cv-foundation.org/openaccess/content_cvpr_2015/papers/Reininghaus_A_Stable_Multi-Scale_2015_CVPR_paper.pdf for more details.
+    :param D1: (n x 2) numpy.array encoding the (finite points of the) first diagram. Must not contain essential points (i.e. with infinite coordinate).
+    :param D2: (m x 2) numpy.array encoding the second diagram.
+    :param bandwidth: bandwidth of the Gaussian kernel with which persistence diagrams will be convolved
+    :param kernel_approx: kernel approximation class used to speed up computation. Common kernel approximations classes can be found in the scikit-learn library (such as RBFSampler for instance).
+    :returns: the persistence scale space kernel value between persistence diagrams.
+    :rtype: float 
+    """
+    DD1 = np.concatenate([D1, D1[:,[1,0]]], axis=0)
+    DD2 = np.concatenate([D2, D2[:,[1,0]]], axis=0)
+    weight_pss = lambda x: 1 if x[1] >= x[0] else -1
+    return 0.5 * persistence_weighted_gaussian_kernel(DD1, DD2, weight=weight_pss, kernel_approx=kernel_approx, bandwidth=bandwidth)
+
+def pairwise_persistence_diagram_kernels(X, Y=None, metric="sliced_wasserstein", **kwargs):
+    """
+    This function computes the kernel matrix between two lists of persistence diagrams given as numpy arrays of shape (nx2).
+    :param X: first list of persistence diagrams. 
+    :param Y: second list of persistence diagrams (optional). If None, pairwise kernel values are computed from the first list only.
+    :param metric: kernel to use. It can be either a string ("sliced_wasserstein", "persistence_scale_space", "persistence_weighted_gaussian", "persistence_fisher") or a function taking two numpy arrays of shape (nx2) and (mx2) as inputs.
+    :returns: kernel matrix, i.e., numpy array of shape (num diagrams 1 x num diagrams 2)
+    :rtype: float
+    """
+    if Y is None:
+        YY = None
+        pX = Padding(use=True).fit_transform(X)
+        diag_len = len(pX[0])
+        XX = np.reshape(np.vstack(pX), [-1, diag_len*3])
+    else:
+        nX, nY = len(X), len(Y)
+        pD = Padding(use=True).fit_transform(X + Y)
+        diag_len = len(pD[0])
+        XX = np.reshape(np.vstack(pD[:nX]), [-1, diag_len*3])
+        YY = np.reshape(np.vstack(pD[nX:]), [-1, diag_len*3])
+
+    if metric == "sliced_wasserstein":
+        return np.exp(-pairwise_persistence_diagram_distances(X, Y, metric="sliced_wasserstein", num_directions=kwargs["num_directions"]) / kwargs["bandwidth"])
+    elif metric == "persistence_fisher":
+        return np.exp(-pairwise_persistence_diagram_distances(X, Y, metric="persistence_fisher", kernel_approx=kwargs["kernel_approx"], bandwidth=kwargs["bandwidth"]) / kwargs["bandwidth_fisher"])
+    elif metric == "persistence_scale_space":
+        return pairwise_kernels(XX, YY, metric=sklearn_wrapper(persistence_scale_space_kernel, **kwargs))
+    elif metric == "persistence_weighted_gaussian":
+        return pairwise_kernels(XX, YY, metric=sklearn_wrapper(persistence_weighted_gaussian_kernel, **kwargs))
+    else:
+        return pairwise_kernels(XX, YY, metric=sklearn_wrapper(metric, **kwargs))
+
 class SlicedWassersteinKernel(BaseEstimator, TransformerMixin):
     """
     This is a class for computing the sliced Wasserstein kernel matrix from a list of persistence diagrams. The sliced Wasserstein kernel is computed by exponentiating the corresponding sliced Wasserstein distance with a Gaussian kernel. See http://proceedings.mlr.press/v70/carriere17a.html for more details. 
@@ -29,7 +99,7 @@ class SlicedWassersteinKernel(BaseEstimator, TransformerMixin):
             num_directions (int): number of lines evenly sampled from [-pi/2,pi/2] in order to approximate and speed up the kernel computation (default 10).
         """
         self.bandwidth = bandwidth
-        self.sw_ = SlicedWassersteinDistance(num_directions=num_directions)
+        self.num_directions = num_directions
 
     def fit(self, X, y=None):
         """
@@ -39,7 +109,7 @@ class SlicedWassersteinKernel(BaseEstimator, TransformerMixin):
             X (list of n x 2 numpy arrays): input persistence diagrams.
             y (n x 1 array): persistence diagram labels (unused).
         """
-        self.sw_.fit(X, y)
+        self.diagrams_ = X
         return self
 
     def transform(self, X):
@@ -52,7 +122,7 @@ class SlicedWassersteinKernel(BaseEstimator, TransformerMixin):
         Returns:
             numpy array of shape (number of diagrams in **diagrams**) x (number of diagrams in X): matrix of pairwise sliced Wasserstein kernel values.
         """
-        return np.exp(-self.sw_.transform(X)/self.bandwidth)
+        return pairwise_persistence_diagram_kernels(X, self.diagrams_, metric="sliced_wasserstein", bandwidth=self.bandwidth, num_directions=self.num_directions) 
 
 class PersistenceWeightedGaussianKernel(BaseEstimator, TransformerMixin):
     """
@@ -78,10 +148,7 @@ class PersistenceWeightedGaussianKernel(BaseEstimator, TransformerMixin):
             X (list of n x 2 numpy arrays): input persistence diagrams.
             y (n x 1 array): persistence diagram labels (unused).
         """
-        self.diagrams_ = list(X)
-        self.ws_ = [ np.array([self.weight(self.diagrams_[i][j,:]) for j in range(self.diagrams_[i].shape[0])]) for i in range(len(self.diagrams_)) ]
-        if self.kernel_approx is not None:
-            self.approx_ = np.concatenate([np.sum(np.multiply(self.ws_[i][:,np.newaxis], self.kernel_approx.transform(self.diagrams_[i])), axis=0)[np.newaxis,:] for i in range(len(self.diagrams_))])
+        self.diagrams_ = X
         return self
 
     def transform(self, X):
@@ -94,31 +161,7 @@ class PersistenceWeightedGaussianKernel(BaseEstimator, TransformerMixin):
         Returns:
             numpy array of shape (number of diagrams in **diagrams**) x (number of diagrams in X): matrix of pairwise persistence weighted Gaussian kernel values.
         """
-        Xp = list(X)
-        Xfit = np.zeros((len(Xp), len(self.diagrams_)))
-        if len(self.diagrams_) == len(Xp) and np.all([np.array_equal(self.diagrams_[i], Xp[i]) for i in range(len(Xp))]):
-            if self.kernel_approx is not None:
-                Xfit = (1./(np.sqrt(2*np.pi)*self.bandwidth)) * np.matmul(self.approx_, self.approx_.T)
-            else:
-                for i in range(len(self.diagrams_)):
-                    for j in range(i+1, len(self.diagrams_)):
-                        W = np.matmul(self.ws_[i][:,np.newaxis], self.ws_[j][np.newaxis,:])
-                        E = (1./(np.sqrt(2*np.pi)*self.bandwidth)) * np.exp(-np.square(pairwise_distances(self.diagrams_[i], self.diagrams_[j]))/(2*np.square(self.bandwidth)))
-                        Xfit[i,j] = np.sum(np.multiply(W, E))
-                        Xfit[j,i] = Xfit[i,j]
-        else:
-            ws = [ np.array([self.weight(Xp[i][j,:]) for j in range(Xp[i].shape[0])]) for i in range(len(Xp)) ]
-            if self.kernel_approx is not None:
-                approx = np.concatenate([np.sum(np.multiply(ws[i][:,np.newaxis], self.kernel_approx.transform(Xp[i])), axis=0)[np.newaxis,:] for i in range(len(Xp))])
-                Xfit = (1./(np.sqrt(2*np.pi)*self.bandwidth)) * np.matmul(approx, self.approx_.T)
-            else:
-                for i in range(len(Xp)):
-                    for j in range(len(self.diagrams_)):
-                        W = np.matmul(ws[i][:,np.newaxis], self.ws_[j][np.newaxis,:])
-                        E = (1./(np.sqrt(2*np.pi)*self.bandwidth)) * np.exp(-np.square(pairwise_distances(Xp[i], self.diagrams_[j]))/(2*np.square(self.bandwidth)))
-                        Xfit[i,j] = np.sum(np.multiply(W, E))
-        
-        return Xfit
+        return pairwise_persistence_diagram_kernels(X, self.diagrams_, metric="persistence_weighted_gaussian", bandwidth=self.bandwidth, weight=self.weight, kernel_approx=self.kernel_approx) 
 
 class PersistenceScaleSpaceKernel(BaseEstimator, TransformerMixin):
     """
@@ -132,7 +175,7 @@ class PersistenceScaleSpaceKernel(BaseEstimator, TransformerMixin):
             bandwidth (double): bandwidth of the Gaussian kernel with which persistence diagrams will be convolved (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.pwg_ = PersistenceWeightedGaussianKernel(bandwidth=bandwidth, weight=lambda x: 1 if x[1] >= x[0] else -1, kernel_approx=kernel_approx)
+        self.bandwidth, self.kernel_approx = bandwidth, kernel_approx
 
     def fit(self, X, y=None):
         """
@@ -142,11 +185,7 @@ class PersistenceScaleSpaceKernel(BaseEstimator, TransformerMixin):
             X (list of n x 2 numpy arrays): input persistence diagrams.
             y (n x 1 array): persistence diagram labels (unused).
         """
-        self.diagrams_ = list(X)
-        for i in range(len(self.diagrams_)):
-            op_D = self.diagrams_[i][:,[1,0]]
-            self.diagrams_[i] = np.concatenate([self.diagrams_[i], op_D], axis=0)
-        self.pwg_.fit(X)
+        self.diagrams_ = X
         return self
 
     def transform(self, X):
@@ -159,11 +198,7 @@ class PersistenceScaleSpaceKernel(BaseEstimator, TransformerMixin):
         Returns:
             numpy array of shape (number of diagrams in **diagrams**) x (number of diagrams in X): matrix of pairwise persistence scale space kernel values.
         """
-        Xp = list(X)
-        for i in range(len(Xp)):
-            op_X = Xp[i][:,[1,0]]
-            Xp[i] = np.concatenate([Xp[i], op_X], axis=0)
-        return self.pwg_.transform(Xp)
+        return pairwise_persistence_diagram_kernels(X, self.diagrams_, metric="persistence_scale_space", bandwidth=self.bandwidth, kernel_approx=self.kernel_approx)
 
 class PersistenceFisherKernel(BaseEstimator, TransformerMixin):
     """
@@ -179,7 +214,7 @@ class PersistenceFisherKernel(BaseEstimator, TransformerMixin):
             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 = bandwidth
-        self.pf_ = PersistenceFisherDistance(bandwidth=bandwidth_fisher, kernel_approx=kernel_approx)
+        self.bandwidth_fisher, self.kernel_approx = bandwidth_fisher, kernel_approx
 
     def fit(self, X, y=None):
         """
@@ -189,7 +224,7 @@ class PersistenceFisherKernel(BaseEstimator, TransformerMixin):
             X (list of n x 2 numpy arrays): input persistence diagrams.
             y (n x 1 array): persistence diagram labels (unused).
         """
-        self.pf_.fit(X, y)
+        self.diagrams_ = X
         return self
 
     def transform(self, X):
@@ -202,5 +237,5 @@ class PersistenceFisherKernel(BaseEstimator, TransformerMixin):
         Returns:
             numpy array of shape (number of diagrams in **diagrams**) x (number of diagrams in X): matrix of pairwise persistence Fisher kernel values.
         """
-        return np.exp(-self.pf_.transform(X)/self.bandwidth)
+        return pairwise_persistence_diagram_kernels(X, self.diagrams_, metric="persistence_fisher", bandwidth=self.bandwidth, bandwidth_fisher=self.bandwidth_fisher, kernel_approx=self.kernel_approx)
 
diff --git a/src/python/gudhi/representations/metrics.py b/src/python/gudhi/representations/metrics.py
index 290c1d07..cc788994 100644
--- a/src/python/gudhi/representations/metrics.py
+++ b/src/python/gudhi/representations/metrics.py
@@ -11,6 +11,8 @@ import numpy as np
 from sklearn.base import BaseEstimator, TransformerMixin
 from sklearn.metrics import pairwise_distances
 from gudhi.wasserstein import wasserstein_distance
+from .preprocessing import Padding
+
 try:
     from .. import bottleneck_distance
     USE_GUDHI = True
@@ -22,6 +24,108 @@ except ImportError:
 # Metrics ###################################
 #############################################
 
+def sliced_wasserstein_distance(D1, D2, num_directions):
+    """
+    This is a function for computing the sliced Wasserstein distance from two 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.
+    :param D1: (n x 2) numpy.array encoding the (finite points of the) first diagram. Must not contain essential points (i.e. with infinite coordinate).
+    :param D2: (m x 2) numpy.array encoding the second diagram.
+    :param num_directions: number of lines evenly sampled from [-pi/2,pi/2] in order to approximate and speed up the distance computation.
+    :returns: the sliced Wasserstein distance between persistence diagrams.
+    :rtype: float 
+    """
+    thetas = np.linspace(-np.pi/2, np.pi/2, num=num_directions+1)[np.newaxis,:-1]
+    lines = np.concatenate([np.cos(thetas), np.sin(thetas)], axis=0)
+    approx1 = np.matmul(D1, lines)
+    diag_proj1 = (1./2) * np.ones((2,2))
+    approx_diag1 = np.matmul(np.matmul(D1, diag_proj1), lines)
+    approx2 = np.matmul(D2, lines)
+    diag_proj2 = (1./2) * np.ones((2,2))
+    approx_diag2 = np.matmul(np.matmul(D2, diag_proj2), lines)
+    A = np.sort(np.concatenate([approx1, approx_diag2], axis=0), axis=0)
+    B = np.sort(np.concatenate([approx2, approx_diag1], axis=0), axis=0)
+    L1 = np.sum(np.abs(A-B), axis=0)
+    return np.mean(L1)
+
+def persistence_fisher_distance(D1, D2, kernel_approx=None, bandwidth=1.):
+    """
+    This is a function for computing the persistence Fisher distance from two 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.
+    :param D1: (n x 2) numpy.array encoding the (finite points of the) first diagram. Must not contain essential points (i.e. with infinite coordinate).
+    :param D2: (m x 2) numpy.array encoding the second diagram.
+    :param bandwidth: bandwidth of the Gaussian kernel used to turn persistence diagrams into probability distributions.
+    :param kernel_approx: kernel approximation class used to speed up computation. Common kernel approximations classes can be found in the scikit-learn library (such as RBFSampler for instance).   
+    :returns: the persistence Fisher distance between persistence diagrams.
+    :rtype: float 
+    """
+    projection = (1./2) * np.ones((2,2))
+    diagonal_projections1 = np.matmul(D1, projection)
+    diagonal_projections2 = np.matmul(D2, projection)
+    if kernel_approx is not None:
+        approx1 = kernel_approx.transform(D1)
+        approx_diagonal1 = kernel_approx.transform(diagonal_projections1)
+        approx2 = kernel_approx.transform(D2)
+        approx_diagonal2 = kernel_approx.transform(diagonal_projections2)
+        Z = np.concatenate([approx1, approx_diagonal1, approx2, approx_diagonal2], axis=0)
+        U, V = np.sum(np.concatenate([approx1, approx_diagonal2], axis=0), axis=0), np.sum(np.concatenate([approx2, approx_diagonal1], 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
+        return np.arccos(  min(np.dot(np.sqrt(vectori), np.sqrt(vectorj)), 1.)  )
+    else:
+        Z = np.concatenate([D1, diagonal_projections1, D2, diagonal_projections2], axis=0)
+        U, V = np.concatenate([D1, diagonal_projections2], axis=0), np.concatenate([D2, diagonal_projections1], axis=0) 
+        vectori = np.sum(np.exp(-np.square(pairwise_distances(Z,U))/(2 * np.square(bandwidth)))/(bandwidth * np.sqrt(2*np.pi)), axis=1)
+        vectorj = np.sum(np.exp(-np.square(pairwise_distances(Z,V))/(2 * np.square(bandwidth)))/(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
+        return np.arccos(  min(np.dot(np.sqrt(vectori), np.sqrt(vectorj)), 1.)  )
+
+def sklearn_wrapper(metric, **kwargs):
+    """
+    This function is a wrapper for any metric between two persistence diagrams that takes two numpy arrays of shapes (nx2) and (mx2) as arguments. It turns the metric into another that takes flattened and padded diagrams as inputs. 
+    """
+    def flat_metric(D1, D2):
+        DD1, DD2 = np.reshape(D1, [-1,3]), np.reshape(D2, [-1,3])
+        return metric(DD1[DD1[:,2]==1,0:2], DD2[DD2[:,2]==1,0:2], **kwargs)
+    return flat_metric
+
+def pairwise_persistence_diagram_distances(X, Y=None, metric="bottleneck", **kwargs):
+    """
+    This function computes the distance matrix between two lists of persistence diagrams given as numpy arrays of shape (nx2).
+    :param X: first list of persistence diagrams. 
+    :param Y: second list of persistence diagrams (optional). If None, pairwise distances are computed from the first list only.
+    :param metric: distance to use. It can be either a string ("sliced_wasserstein", "wasserstein", "bottleneck", "persistence_fisher") or a function taking two numpy arrays of shape (nx2) and (mx2) as inputs.
+    :returns: distance matrix, i.e., numpy array of shape (num diagrams 1 x num diagrams 2)
+    :rtype: float
+    """
+    if Y is None:
+        YY = None
+        pX = Padding(use=True).fit_transform(X)
+        diag_len = len(pX[0])
+        XX = np.reshape(np.vstack(pX), [-1, diag_len*3])
+    else:
+        nX, nY = len(X), len(Y)
+        pD = Padding(use=True).fit_transform(X + Y)
+        diag_len = len(pD[0])
+        XX = np.reshape(np.vstack(pD[:nX]), [-1, diag_len*3])
+        YY = np.reshape(np.vstack(pD[nX:]), [-1, diag_len*3])
+
+    if metric == "bottleneck":
+        return pairwise_distances(XX, YY, metric=sklearn_wrapper(bottleneck_distance, **kwargs))
+    elif metric == "wasserstein":
+        return pairwise_distances(XX, YY, metric=sklearn_wrapper(wasserstein_distance, **kwargs))
+    elif metric == "sliced_wasserstein":
+        return pairwise_distances(XX, YY, metric=sklearn_wrapper(sliced_wasserstein_distance, **kwargs))
+    elif metric == "persistence_fisher":
+        return pairwise_distances(XX, YY, metric=sklearn_wrapper(persistence_fisher_distance, **kwargs))
+    else:
+        return pairwise_distances(XX, YY, metric=sklearn_wrapper(metric, **kwargs))
+
 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. 
@@ -34,8 +138,6 @@ class SlicedWassersteinDistance(BaseEstimator, TransformerMixin):
             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):
         """
@@ -46,9 +148,6 @@ class SlicedWassersteinDistance(BaseEstimator, TransformerMixin):
             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):
@@ -61,27 +160,7 @@ class SlicedWassersteinDistance(BaseEstimator, TransformerMixin):
         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
+        return pairwise_persistence_diagram_distances(X, self.diagrams_, metric="sliced_wasserstein", num_directions=self.num_directions)
 
 class BottleneckDistance(BaseEstimator, TransformerMixin):
     """
@@ -117,33 +196,9 @@ class BottleneckDistance(BaseEstimator, TransformerMixin):
         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
-
+        if not USE_GUDHI:
+            print("Gudhi built without CGAL: returning a null matrix")
+        Xfit = pairwise_persistence_diagram_distances(X, self.diagrams_, metric="bottleneck", e=self.epsilon) if USE_GUDHI else np.zeros((len(X), len(self.diagrams_)))
         return Xfit
 
 class WassersteinDistance(BaseEstimator, TransformerMixin):
@@ -181,28 +236,7 @@ class WassersteinDistance(BaseEstimator, TransformerMixin):
         Returns:
             numpy array of shape (number of diagrams in **diagrams**) x (number of diagrams in X): matrix of pairwise Wasserstein 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))
-
-            for i in range(num_diag1):
-                for j in range(i+1, num_diag1):
-                    matrix[i,j] = wasserstein_distance(X[i], X[j], self.order, self.internal_p)
-                    matrix[j,i] = matrix[i,j]
-
-        else:
-            num_diag2 = len(self.diagrams_)
-            matrix = np.zeros((num_diag1, num_diag2))
-
-            for i in range(num_diag1):
-                for j in range(num_diag2):
-                    matrix[i,j] = wasserstein_distance(X[i], self.diagrams_[j], self.order, self.internal_p)
-
-        Xfit = matrix 
-
-        return Xfit
+        return pairwise_persistence_diagram_distances(X, self.diagrams_, metric="wasserstein", order=self.order, internal_p=self.internal_p)
 
 class PersistenceFisherDistance(BaseEstimator, TransformerMixin):
     """
@@ -227,11 +261,6 @@ class PersistenceFisherDistance(BaseEstimator, TransformerMixin):
             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):
@@ -244,60 +273,4 @@ class PersistenceFisherDistance(BaseEstimator, TransformerMixin):
         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
+        return pairwise_persistence_diagram_distances(X, self.diagrams_, metric="persistence_fisher", bandwidth=self.bandwidth, kernel_approx=self.kernel_approx)
-- 
cgit v1.2.3


From 29e81d5038116aef0ec505e4d21d29f1c5920e34 Mon Sep 17 00:00:00 2001
From: MathieuCarriere <mathieu.carriere3@gmail.com>
Date: Fri, 7 Feb 2020 21:00:17 -0500
Subject: added sklearn trick

---
 src/python/gudhi/representations/kernel_methods.py | 20 +++---------
 src/python/gudhi/representations/metrics.py        | 37 +++++++++-------------
 2 files changed, 20 insertions(+), 37 deletions(-)

diff --git a/src/python/gudhi/representations/kernel_methods.py b/src/python/gudhi/representations/kernel_methods.py
index bbbb7c31..d89f69ab 100644
--- a/src/python/gudhi/representations/kernel_methods.py
+++ b/src/python/gudhi/representations/kernel_methods.py
@@ -62,27 +62,17 @@ def pairwise_persistence_diagram_kernels(X, Y=None, metric="sliced_wasserstein",
     :param metric: kernel to use. It can be either a string ("sliced_wasserstein", "persistence_scale_space", "persistence_weighted_gaussian", "persistence_fisher") or a function taking two numpy arrays of shape (nx2) and (mx2) as inputs.
     :returns: kernel matrix, i.e., numpy array of shape (num diagrams 1 x num diagrams 2)
     :rtype: float
-    """
-    if Y is None:
-        YY = None
-        pX = Padding(use=True).fit_transform(X)
-        diag_len = len(pX[0])
-        XX = np.reshape(np.vstack(pX), [-1, diag_len*3])
-    else:
-        nX, nY = len(X), len(Y)
-        pD = Padding(use=True).fit_transform(X + Y)
-        diag_len = len(pD[0])
-        XX = np.reshape(np.vstack(pD[:nX]), [-1, diag_len*3])
-        YY = np.reshape(np.vstack(pD[nX:]), [-1, diag_len*3])
-
+    """    
+    XX = np.reshape(np.arange(len(X)), [-1,1])
+    YY = None if Y is None else np.reshape(np.arange(len(Y)), [-1,1])
     if metric == "sliced_wasserstein":
         return np.exp(-pairwise_persistence_diagram_distances(X, Y, metric="sliced_wasserstein", num_directions=kwargs["num_directions"]) / kwargs["bandwidth"])
     elif metric == "persistence_fisher":
         return np.exp(-pairwise_persistence_diagram_distances(X, Y, metric="persistence_fisher", kernel_approx=kwargs["kernel_approx"], bandwidth=kwargs["bandwidth"]) / kwargs["bandwidth_fisher"])
     elif metric == "persistence_scale_space":
-        return pairwise_kernels(XX, YY, metric=sklearn_wrapper(persistence_scale_space_kernel, **kwargs))
+        return pairwise_kernels(XX, YY, metric=sklearn_wrapper(persistence_scale_space_kernel, X, Y, **kwargs))
     elif metric == "persistence_weighted_gaussian":
-        return pairwise_kernels(XX, YY, metric=sklearn_wrapper(persistence_weighted_gaussian_kernel, **kwargs))
+        return pairwise_kernels(XX, YY, metric=sklearn_wrapper(persistence_weighted_gaussian_kernel, X, Y, **kwargs))
     else:
         return pairwise_kernels(XX, YY, metric=sklearn_wrapper(metric, **kwargs))
 
diff --git a/src/python/gudhi/representations/metrics.py b/src/python/gudhi/representations/metrics.py
index cc788994..fead8aa0 100644
--- a/src/python/gudhi/representations/metrics.py
+++ b/src/python/gudhi/representations/metrics.py
@@ -85,13 +85,16 @@ def persistence_fisher_distance(D1, D2, kernel_approx=None, bandwidth=1.):
             vectorj = vectorj/vectorj_sum
         return np.arccos(  min(np.dot(np.sqrt(vectori), np.sqrt(vectorj)), 1.)  )
 
-def sklearn_wrapper(metric, **kwargs):
+def sklearn_wrapper(metric, X, Y, **kwargs):
     """
-    This function is a wrapper for any metric between two persistence diagrams that takes two numpy arrays of shapes (nx2) and (mx2) as arguments. It turns the metric into another that takes flattened and padded diagrams as inputs. 
+    This function is a wrapper for any metric between two persistence diagrams that takes two numpy arrays of shapes (nx2) and (mx2) as arguments.
     """
-    def flat_metric(D1, D2):
-        DD1, DD2 = np.reshape(D1, [-1,3]), np.reshape(D2, [-1,3])
-        return metric(DD1[DD1[:,2]==1,0:2], DD2[DD2[:,2]==1,0:2], **kwargs)
+    if Y is None:
+        def flat_metric(a, b):
+            return metric(X[int(a[0])], X[int(b[0])], **kwargs)
+    else:
+        def flat_metric(a, b):
+            return metric(X[int(a[0])], Y[int(b[0])], **kwargs)
     return flat_metric
 
 def pairwise_persistence_diagram_distances(X, Y=None, metric="bottleneck", **kwargs):
@@ -103,28 +106,18 @@ def pairwise_persistence_diagram_distances(X, Y=None, metric="bottleneck", **kwa
     :returns: distance matrix, i.e., numpy array of shape (num diagrams 1 x num diagrams 2)
     :rtype: float
     """
-    if Y is None:
-        YY = None
-        pX = Padding(use=True).fit_transform(X)
-        diag_len = len(pX[0])
-        XX = np.reshape(np.vstack(pX), [-1, diag_len*3])
-    else:
-        nX, nY = len(X), len(Y)
-        pD = Padding(use=True).fit_transform(X + Y)
-        diag_len = len(pD[0])
-        XX = np.reshape(np.vstack(pD[:nX]), [-1, diag_len*3])
-        YY = np.reshape(np.vstack(pD[nX:]), [-1, diag_len*3])
-
+    XX = np.reshape(np.arange(len(X)), [-1,1])
+    YY = None if Y is None else np.reshape(np.arange(len(Y)), [-1,1]) 
     if metric == "bottleneck":
-        return pairwise_distances(XX, YY, metric=sklearn_wrapper(bottleneck_distance, **kwargs))
+        return pairwise_distances(XX, YY, metric=sklearn_wrapper(bottleneck_distance, X, Y, **kwargs))
     elif metric == "wasserstein":
-        return pairwise_distances(XX, YY, metric=sklearn_wrapper(wasserstein_distance, **kwargs))
+        return pairwise_distances(XX, YY, metric=sklearn_wrapper(wasserstein_distance, X, Y, **kwargs))
     elif metric == "sliced_wasserstein":
-        return pairwise_distances(XX, YY, metric=sklearn_wrapper(sliced_wasserstein_distance, **kwargs))
+        return pairwise_distances(XX, YY, metric=sklearn_wrapper(sliced_wasserstein_distance, X, Y, **kwargs))
     elif metric == "persistence_fisher":
-        return pairwise_distances(XX, YY, metric=sklearn_wrapper(persistence_fisher_distance, **kwargs))
+        return pairwise_distances(XX, YY, metric=sklearn_wrapper(persistence_fisher_distance, X, Y, **kwargs))
     else:
-        return pairwise_distances(XX, YY, metric=sklearn_wrapper(metric, **kwargs))
+        return pairwise_distances(XX, YY, metric=sklearn_wrapper(metric, X, Y, **kwargs))
 
 class SlicedWassersteinDistance(BaseEstimator, TransformerMixin):
     """
-- 
cgit v1.2.3


From ef0f82ef2155440827e17c552abb49b509866fc7 Mon Sep 17 00:00:00 2001
From: mathieu <mathieu.carriere3@gmail.com>
Date: Thu, 13 Feb 2020 16:01:29 -0500
Subject: integrated hera

---
 .../diagram_vectorizations_distances_kernels.py    |  7 ++++++-
 src/python/gudhi/representations/metrics.py        | 23 ++++++++++++++++------
 2 files changed, 23 insertions(+), 7 deletions(-)

diff --git a/src/python/example/diagram_vectorizations_distances_kernels.py b/src/python/example/diagram_vectorizations_distances_kernels.py
index 66c32cc2..6352d2b5 100755
--- a/src/python/example/diagram_vectorizations_distances_kernels.py
+++ b/src/python/example/diagram_vectorizations_distances_kernels.py
@@ -117,7 +117,12 @@ X = SW.fit(diags)
 Y = SW.transform(diags2)
 print("SW kernel is " + str(Y[0][0]))
 
-W = WassersteinDistance(order=2, internal_p=2)
+W = WassersteinDistance(order=2, internal_p=2, mode="pot")
+X = W.fit(diags)
+Y = W.transform(diags2)
+print("Wasserstein distance is " + str(Y[0][0]))
+
+W = WassersteinDistance(order=2, internal_p=2, mode="hera", delta=0.0001)
 X = W.fit(diags)
 Y = W.transform(diags2)
 print("Wasserstein distance is " + str(Y[0][0]))
diff --git a/src/python/gudhi/representations/metrics.py b/src/python/gudhi/representations/metrics.py
index cc788994..ed998603 100644
--- a/src/python/gudhi/representations/metrics.py
+++ b/src/python/gudhi/representations/metrics.py
@@ -10,7 +10,8 @@
 import numpy as np
 from sklearn.base import BaseEstimator, TransformerMixin
 from sklearn.metrics import pairwise_distances
-from gudhi.wasserstein import wasserstein_distance
+from gudhi.wasserstein import wasserstein_distance as pot_wasserstein_distance
+from gudhi.hera import wasserstein_distance as hera_wasserstein_distance
 from .preprocessing import Padding
 
 try:
@@ -117,8 +118,10 @@ def pairwise_persistence_diagram_distances(X, Y=None, metric="bottleneck", **kwa
 
     if metric == "bottleneck":
         return pairwise_distances(XX, YY, metric=sklearn_wrapper(bottleneck_distance, **kwargs))
-    elif metric == "wasserstein":
-        return pairwise_distances(XX, YY, metric=sklearn_wrapper(wasserstein_distance, **kwargs))
+    elif metric == "wasserstein" or metric == "pot_wasserstein":
+        return pairwise_distances(XX, YY, metric=sklearn_wrapper(pot_wasserstein_distance, **kwargs))
+    elif metric == "hera_wasserstein":
+        return pairwise_distances(XX, YY, metric=sklearn_wrapper(hera_wasserstein_distance, **kwargs))
     elif metric == "sliced_wasserstein":
         return pairwise_distances(XX, YY, metric=sklearn_wrapper(sliced_wasserstein_distance, **kwargs))
     elif metric == "persistence_fisher":
@@ -205,15 +208,19 @@ class WassersteinDistance(BaseEstimator, TransformerMixin):
     """
     This is a class for computing the Wasserstein distance matrix from a list of persistence diagrams. 
     """
-    def __init__(self, order=2, internal_p=2):
+    def __init__(self, order=2, internal_p=2, mode="pot", delta=0.0001):
         """
         Constructor for the WassersteinDistance class.
 
         Parameters:
             order (int): exponent for Wasserstein, default value is 2., see :func:`gudhi.wasserstein.wasserstein_distance`.
             internal_p (int): ground metric on the (upper-half) plane (i.e. norm l_p in R^2), default value is 2 (euclidean norm), see :func:`gudhi.wasserstein.wasserstein_distance`.
+            mode (str): method for computing Wasserstein distance. Either "pot" or "hera".
+            delta (float): relative error 1+delta. Used only if mode == "hera".
         """
-        self.order, self.internal_p = order, internal_p
+        self.order, self.internal_p, self.mode = order, internal_p, mode
+        self.metric = "pot_wasserstein" if mode == "pot" else "hera_wasserstein"
+        self.delta = delta
 
     def fit(self, X, y=None):
         """
@@ -236,7 +243,11 @@ class WassersteinDistance(BaseEstimator, TransformerMixin):
         Returns:
             numpy array of shape (number of diagrams in **diagrams**) x (number of diagrams in X): matrix of pairwise Wasserstein distances.
         """
-        return pairwise_persistence_diagram_distances(X, self.diagrams_, metric="wasserstein", order=self.order, internal_p=self.internal_p)
+        if self.metric == "hera_wasserstein":
+            Xfit = pairwise_persistence_diagram_distances(X, self.diagrams_, metric=self.metric, order=self.order, internal_p=self.internal_p, delta=self.delta)
+        else:
+            Xfit = pairwise_persistence_diagram_distances(X, self.diagrams_, metric=self.metric, order=self.order, internal_p=self.internal_p)
+        return Xfit
 
 class PersistenceFisherDistance(BaseEstimator, TransformerMixin):
     """
-- 
cgit v1.2.3


From d9290a78741fc14dc0f87d395da967a4d561b34a Mon Sep 17 00:00:00 2001
From: mathieu <mathieu.carriere3@gmail.com>
Date: Thu, 13 Feb 2020 16:11:34 -0500
Subject: small modif on example file

---
 src/python/example/diagram_vectorizations_distances_kernels.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/python/example/diagram_vectorizations_distances_kernels.py b/src/python/example/diagram_vectorizations_distances_kernels.py
index 6352d2b5..507ead7c 100755
--- a/src/python/example/diagram_vectorizations_distances_kernels.py
+++ b/src/python/example/diagram_vectorizations_distances_kernels.py
@@ -120,12 +120,12 @@ print("SW kernel is " + str(Y[0][0]))
 W = WassersteinDistance(order=2, internal_p=2, mode="pot")
 X = W.fit(diags)
 Y = W.transform(diags2)
-print("Wasserstein distance is " + str(Y[0][0]))
+print("Wasserstein distance (POT) is " + str(Y[0][0]))
 
 W = WassersteinDistance(order=2, internal_p=2, mode="hera", delta=0.0001)
 X = W.fit(diags)
 Y = W.transform(diags2)
-print("Wasserstein distance is " + str(Y[0][0]))
+print("Wasserstein distance (hera) is " + str(Y[0][0]))
 
 W = BottleneckDistance(epsilon=.001)
 X = W.fit(diags)
-- 
cgit v1.2.3


From a47ace987876cb52351ae9223d335629aedbd71e Mon Sep 17 00:00:00 2001
From: mathieu <mathieu.carriere3@gmail.com>
Date: Tue, 10 Mar 2020 19:44:57 -0400
Subject: new fixes

---
 ext/hera                                    |  2 +-
 src/python/gudhi/representations/metrics.py | 27 ++++++++++++---------------
 2 files changed, 13 insertions(+), 16 deletions(-)

diff --git a/ext/hera b/ext/hera
index cb1838e6..9a899718 160000
--- a/ext/hera
+++ b/ext/hera
@@ -1 +1 @@
-Subproject commit cb1838e682ec07f80720241cf9098400caeb83c7
+Subproject commit 9a89971855acefe39dce0e2adadf53b88ca8f683
diff --git a/src/python/gudhi/representations/metrics.py b/src/python/gudhi/representations/metrics.py
index c5439a67..0659b457 100644
--- a/src/python/gudhi/representations/metrics.py
+++ b/src/python/gudhi/representations/metrics.py
@@ -10,17 +10,9 @@
 import numpy as np
 from sklearn.base import BaseEstimator, TransformerMixin
 from sklearn.metrics import pairwise_distances
-from gudhi.wasserstein import wasserstein_distance as pot_wasserstein_distance
 from gudhi.hera import wasserstein_distance as hera_wasserstein_distance
 from .preprocessing import Padding
 
-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 ###################################
 #############################################
@@ -111,9 +103,13 @@ def pairwise_persistence_diagram_distances(X, Y=None, metric="bottleneck", **kwa
     YY = None if Y is None else np.reshape(np.arange(len(Y)), [-1,1]) 
     if metric == "bottleneck":
         return pairwise_distances(XX, YY, metric=sklearn_wrapper(bottleneck_distance, X, Y, **kwargs))
-    elif metric == "wasserstein" or metric == "pot_wasserstein":
-        return pairwise_distances(XX, YY, metric=sklearn_wrapper(pot_wasserstein_distance,  X, Y, **kwargs))
-    elif metric == "hera_wasserstein":
+    elif metric == "pot_wasserstein":
+        try:
+            from gudhi.wasserstein import wasserstein_distance as pot_wasserstein_distance
+            return pairwise_distances(XX, YY, metric=sklearn_wrapper(pot_wasserstein_distance,  X, Y, **kwargs))
+        except ImportError:
+            print("Gudhi built without POT")
+    elif metric == "wasserstein" or metric == "hera_wasserstein":
         return pairwise_distances(XX, YY, metric=sklearn_wrapper(hera_wasserstein_distance,  X, Y, **kwargs))
     elif metric == "sliced_wasserstein":
         return pairwise_distances(XX, YY, metric=sklearn_wrapper(sliced_wasserstein_distance, X, Y, **kwargs))
@@ -192,16 +188,17 @@ class BottleneckDistance(BaseEstimator, TransformerMixin):
         Returns:
             numpy array of shape (number of diagrams in **diagrams**) x (number of diagrams in X): matrix of pairwise bottleneck distances.
         """
-        if not USE_GUDHI:
-            print("Gudhi built without CGAL: returning a null matrix")
-        Xfit = pairwise_persistence_diagram_distances(X, self.diagrams_, metric="bottleneck", e=self.epsilon) if USE_GUDHI else np.zeros((len(X), len(self.diagrams_)))
+        try:
+            Xfit = pairwise_persistence_diagram_distances(X, self.diagrams_, metric="bottleneck", e=self.epsilon)
+        except ImportError:
+            print("Gudhi built without CGAL")
         return Xfit
 
 class WassersteinDistance(BaseEstimator, TransformerMixin):
     """
     This is a class for computing the Wasserstein distance matrix from a list of persistence diagrams. 
     """
-    def __init__(self, order=2, internal_p=2, mode="pot", delta=0.0001):
+    def __init__(self, order=2, internal_p=2, mode="pot", delta=0.01):
         """
         Constructor for the WassersteinDistance class.
 
-- 
cgit v1.2.3


From 25e40a52ec7bc9e1bfe418fb1aa16e2a06994d1b Mon Sep 17 00:00:00 2001
From: mathieu <mathieu.carriere3@gmail.com>
Date: Wed, 11 Mar 2020 15:35:37 -0400
Subject: new fixes

---
 src/python/gudhi/representations/metrics.py | 63 +++++++++++++++++++++++------
 1 file changed, 50 insertions(+), 13 deletions(-)

diff --git a/src/python/gudhi/representations/metrics.py b/src/python/gudhi/representations/metrics.py
index 0659b457..f913f1fc 100644
--- a/src/python/gudhi/representations/metrics.py
+++ b/src/python/gudhi/representations/metrics.py
@@ -19,7 +19,7 @@ from .preprocessing import Padding
 
 def sliced_wasserstein_distance(D1, D2, num_directions):
     """
-    This is a function for computing the sliced Wasserstein distance from two 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.
+    This is a function for computing the sliced Wasserstein distance from two persistence diagrams. The Sliced Wasserstein distance is computed by projecting the persistence diagrams onto lines, comparing the projections with the 1-norm, and finally averaging over the lines. See http://proceedings.mlr.press/v70/carriere17a.html for more details.
     :param D1: (n x 2) numpy.array encoding the (finite points of the) first diagram. Must not contain essential points (i.e. with infinite coordinate).
     :param D2: (m x 2) numpy.array encoding the second diagram.
     :param num_directions: number of lines evenly sampled from [-pi/2,pi/2] in order to approximate and speed up the distance computation.
@@ -39,6 +39,34 @@ def sliced_wasserstein_distance(D1, D2, num_directions):
     L1 = np.sum(np.abs(A-B), axis=0)
     return np.mean(L1)
 
+def compute_persistence_diagram_projections(X, num_directions):
+    """
+    This is a function for projecting the points of a list of persistence diagrams (as well as their diagonal projections) onto a fixed number of lines sampled uniformly on [-pi/2, pi/2]. This function can be used as a preprocessing step in order to speed up the running time for computing all pairwise sliced Wasserstein distances / kernel values on a list of persistence diagrams. 
+    :param X: list of persistence diagrams. 
+    :param num_directions: number of lines evenly sampled from [-pi/2,pi/2] in order to approximate and speed up the distance computation.
+    :returns: list of projected persistence diagrams.
+    :rtype: float
+    """
+    thetas = np.linspace(-np.pi/2, np.pi/2, num=num_directions+1)[np.newaxis,:-1]
+    lines = np.concatenate([np.cos(thetas), np.sin(thetas)], axis=0)
+    XX = [np.vstack([np.matmul(D, lines), np.matmul(np.matmul(D, .5 * np.ones((2,2))), lines)]) for D in X]
+    return XX
+
+def sliced_wasserstein_distance_on_projections(D1, D2):
+    """
+    This is a function for computing the sliced Wasserstein distance between two persistence diagrams that have already been projected onto some lines. It simply amounts to comparing the sorted projections with the 1-norm, and averaging over the lines. See http://proceedings.mlr.press/v70/carriere17a.html for more details.
+    :param D1: (2n x number_of_lines) numpy.array containing the n projected points of the first diagram, and the n projections of their diagonal projections.
+    :param D2: (2m x number_of_lines) numpy.array containing the m projected points of the second diagram, and the m projections of their diagonal projections.
+    :returns: the sliced Wasserstein distance between the projected persistence diagrams.
+    :rtype: float 
+    """
+    lim1, lim2 = int(len(D1)/2), int(len(D2)/2)
+    approx1, approx_diag1, approx2, approx_diag2 = D1[:lim1], D1[lim1:], D2[:lim2], D2[lim2:]
+    A = np.sort(np.concatenate([approx1, approx_diag2], axis=0), axis=0)
+    B = np.sort(np.concatenate([approx2, approx_diag1], axis=0), axis=0)
+    L1 = np.sum(np.abs(A-B), axis=0)
+    return np.mean(L1)
+
 def persistence_fisher_distance(D1, D2, kernel_approx=None, bandwidth=1.):
     """
     This is a function for computing the persistence Fisher distance from two 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.
@@ -90,31 +118,43 @@ def sklearn_wrapper(metric, X, Y, **kwargs):
             return metric(X[int(a[0])], Y[int(b[0])], **kwargs)
     return flat_metric
 
+PAIRWISE_DISTANCE_FUNCTIONS = {
+    "wasserstein": hera_wasserstein_distance,
+    "hera_wasserstein": hera_wasserstein_distance,
+    "persistence_fisher": persistence_fisher_distance,
+}
+
 def pairwise_persistence_diagram_distances(X, Y=None, metric="bottleneck", **kwargs):
     """
     This function computes the distance matrix between two lists of persistence diagrams given as numpy arrays of shape (nx2).
     :param X: first list of persistence diagrams. 
     :param Y: second list of persistence diagrams (optional). If None, pairwise distances are computed from the first list only.
-    :param metric: distance to use. It can be either a string ("sliced_wasserstein", "wasserstein", "bottleneck", "persistence_fisher") or a function taking two numpy arrays of shape (nx2) and (mx2) as inputs.
+    :param metric: distance to use. It can be either a string ("sliced_wasserstein", "wasserstein", "hera_wasserstein" (Wasserstein distance computed with Hera---note that Hera is also used for the default option "wasserstein"), "pot_wasserstein" (Wasserstein distance computed with POT), "bottleneck", "persistence_fisher") or a function taking two numpy arrays of shape (nx2) and (mx2) as inputs.
     :returns: distance matrix, i.e., numpy array of shape (num diagrams 1 x num diagrams 2)
     :rtype: float
     """
     XX = np.reshape(np.arange(len(X)), [-1,1])
     YY = None if Y is None else np.reshape(np.arange(len(Y)), [-1,1]) 
     if metric == "bottleneck":
-        return pairwise_distances(XX, YY, metric=sklearn_wrapper(bottleneck_distance, X, Y, **kwargs))
+        try: 
+            from .. import bottleneck_distance
+            return pairwise_distances(XX, YY, metric=sklearn_wrapper(bottleneck_distance, X, Y, **kwargs))
+        except ImportError:
+            print("Gudhi built without CGAL")
+            raise
     elif metric == "pot_wasserstein":
         try:
             from gudhi.wasserstein import wasserstein_distance as pot_wasserstein_distance
             return pairwise_distances(XX, YY, metric=sklearn_wrapper(pot_wasserstein_distance,  X, Y, **kwargs))
         except ImportError:
-            print("Gudhi built without POT")
-    elif metric == "wasserstein" or metric == "hera_wasserstein":
-        return pairwise_distances(XX, YY, metric=sklearn_wrapper(hera_wasserstein_distance,  X, Y, **kwargs))
+            print("Gudhi built without POT. Please install POT or use metric='wasserstein' or metric='hera_wasserstein'")
+            raise
     elif metric == "sliced_wasserstein":
-        return pairwise_distances(XX, YY, metric=sklearn_wrapper(sliced_wasserstein_distance, X, Y, **kwargs))
-    elif metric == "persistence_fisher":
-        return pairwise_distances(XX, YY, metric=sklearn_wrapper(persistence_fisher_distance, X, Y, **kwargs))
+        Xproj = compute_persistence_diagram_projections(X, **kwargs)
+        Yproj = None if Y is None else compute_persistence_diagram_projections(Y, **kwargs)
+        return pairwise_distances(XX, YY, metric=sklearn_wrapper(sliced_wasserstein_distance_on_projections, Xproj, Yproj))
+    elif type(metric) == str:
+        return pairwise_distances(XX, YY, metric=sklearn_wrapper(PAIRWISE_DISTANCE_FUNCTIONS[metric], X, Y, **kwargs))
     else:
         return pairwise_distances(XX, YY, metric=sklearn_wrapper(metric, X, Y, **kwargs))
 
@@ -188,10 +228,7 @@ class BottleneckDistance(BaseEstimator, TransformerMixin):
         Returns:
             numpy array of shape (number of diagrams in **diagrams**) x (number of diagrams in X): matrix of pairwise bottleneck distances.
         """
-        try:
-            Xfit = pairwise_persistence_diagram_distances(X, self.diagrams_, metric="bottleneck", e=self.epsilon)
-        except ImportError:
-            print("Gudhi built without CGAL")
+        Xfit = pairwise_persistence_diagram_distances(X, self.diagrams_, metric="bottleneck", e=self.epsilon)
         return Xfit
 
 class WassersteinDistance(BaseEstimator, TransformerMixin):
-- 
cgit v1.2.3


From 6552d09c3f290a25ee910e007084fe3809f8c8ed Mon Sep 17 00:00:00 2001
From: MathieuCarriere <mathieu.carriere3@gmail.com>
Date: Thu, 12 Mar 2020 16:19:34 -0400
Subject: fixed error message

---
 src/python/gudhi/representations/metrics.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/python/gudhi/representations/metrics.py b/src/python/gudhi/representations/metrics.py
index f913f1fc..4070c321 100644
--- a/src/python/gudhi/representations/metrics.py
+++ b/src/python/gudhi/representations/metrics.py
@@ -147,7 +147,7 @@ def pairwise_persistence_diagram_distances(X, Y=None, metric="bottleneck", **kwa
             from gudhi.wasserstein import wasserstein_distance as pot_wasserstein_distance
             return pairwise_distances(XX, YY, metric=sklearn_wrapper(pot_wasserstein_distance,  X, Y, **kwargs))
         except ImportError:
-            print("Gudhi built without POT. Please install POT or use metric='wasserstein' or metric='hera_wasserstein'")
+            print("POT (Python Optimal Transport) is not installed. Please install POT or use metric='wasserstein' or metric='hera_wasserstein'")
             raise
     elif metric == "sliced_wasserstein":
         Xproj = compute_persistence_diagram_projections(X, **kwargs)
-- 
cgit v1.2.3


From d239af744539572b485b09031f60121383fc1bc6 Mon Sep 17 00:00:00 2001
From: mathieu <mathieu.carriere3@gmail.com>
Date: Fri, 13 Mar 2020 11:31:19 -0400
Subject: tried to fix hera's update

---
 ext/hera | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/ext/hera b/ext/hera
index 9a899718..cb1838e6 160000
--- a/ext/hera
+++ b/ext/hera
@@ -1 +1 @@
-Subproject commit 9a89971855acefe39dce0e2adadf53b88ca8f683
+Subproject commit cb1838e682ec07f80720241cf9098400caeb83c7
-- 
cgit v1.2.3


From 4fb8ab586088dd582b3949cecc11395c37b6f3e6 Mon Sep 17 00:00:00 2001
From: mathieu <mathieu.carriere3@gmail.com>
Date: Fri, 13 Mar 2020 11:32:22 -0400
Subject: tried to fix hera's update

---
 ext/hera | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/ext/hera b/ext/hera
index 0019cae9..cb1838e6 160000
--- a/ext/hera
+++ b/ext/hera
@@ -1 +1 @@
-Subproject commit 0019cae9dc1e9d11aa03bc59681435ba7f21eea8
+Subproject commit cb1838e682ec07f80720241cf9098400caeb83c7
-- 
cgit v1.2.3


From 6410abe3788e17a24b1569bcd7f121d126e1c6cc Mon Sep 17 00:00:00 2001
From: mathieu <mathieu.carriere3@gmail.com>
Date: Fri, 13 Mar 2020 11:58:25 -0400
Subject: fix hera

---
 ext/hera | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/ext/hera b/ext/hera
index cb1838e6..0019cae9 160000
--- a/ext/hera
+++ b/ext/hera
@@ -1 +1 @@
-Subproject commit cb1838e682ec07f80720241cf9098400caeb83c7
+Subproject commit 0019cae9dc1e9d11aa03bc59681435ba7f21eea8
-- 
cgit v1.2.3


From 3099b2395fa143aa6c9b3df2c6087ccd017ff87c Mon Sep 17 00:00:00 2001
From: MathieuCarriere <mathieu.carriere3@gmail.com>
Date: Mon, 16 Mar 2020 12:51:34 -0400
Subject: fixed doc

---
 src/python/gudhi/representations/kernel_methods.py | 45 +++++++++-------
 src/python/gudhi/representations/metrics.py        | 63 +++++++++++++---------
 2 files changed, 66 insertions(+), 42 deletions(-)

diff --git a/src/python/gudhi/representations/kernel_methods.py b/src/python/gudhi/representations/kernel_methods.py
index d89f69ab..50186d63 100644
--- a/src/python/gudhi/representations/kernel_methods.py
+++ b/src/python/gudhi/representations/kernel_methods.py
@@ -20,13 +20,16 @@ from .preprocessing import Padding
 def persistence_weighted_gaussian_kernel(D1, D2, weight=lambda x: 1, kernel_approx=None, bandwidth=1.):
     """
     This is a function for computing the persistence weighted Gaussian kernel value from two persistence diagrams. The persistence weighted Gaussian kernel is computed by convolving the persistence diagram points with weighted Gaussian kernels. See http://proceedings.mlr.press/v48/kusano16.html for more details.
-    :param D1: (n x 2) numpy.array encoding the (finite points of the) first diagram. Must not contain essential points (i.e. with infinite coordinate).
-    :param D2: (m x 2) numpy.array encoding the second diagram.
-    :param bandwidth: bandwidth of the Gaussian kernel with which persistence diagrams will be convolved
-    :param weight: weight function for the persistence diagram points. This function must be defined on 2D points, ie lists or numpy arrays of the form [p_x,p_y].
-    :param kernel_approx: kernel approximation class used to speed up computation. Common kernel approximations classes can be found in the scikit-learn library (such as RBFSampler for instance).
-    :returns: the persistence weighted Gaussian kernel value between persistence diagrams.
-    :rtype: float 
+
+    Parameters:
+        D1: (n x 2) numpy.array encoding the (finite points of the) first diagram. Must not contain essential points (i.e. with infinite coordinate).
+        D2: (m x 2) numpy.array encoding the second diagram.
+        bandwidth (double): bandwidth of the Gaussian kernel with which persistence diagrams will be convolved
+        weight: weight function for the persistence diagram points. This function must be defined on 2D points, ie lists or numpy arrays of the form [p_x,p_y].
+        kernel_approx: kernel approximation class used to speed up computation. Common kernel approximations classes can be found in the scikit-learn library (such as RBFSampler for instance).
+
+    Returns:
+        float: the persistence weighted Gaussian kernel value between persistence diagrams. 
     """
     ws1 = np.array([weight(D1[j,:]) for j in range(len(D1))])
     ws2 = np.array([weight(D2[j,:]) for j in range(len(D2))])
@@ -42,12 +45,15 @@ def persistence_weighted_gaussian_kernel(D1, D2, weight=lambda x: 1, kernel_appr
 def persistence_scale_space_kernel(D1, D2, kernel_approx=None, bandwidth=1.):
     """
     This is a function for computing the persistence scale space kernel value from two persistence diagrams. The persistence scale space kernel is computed by adding the symmetric to the diagonal of each point in each persistence diagram, with negative weight, and then convolving the points with a Gaussian kernel. See https://www.cv-foundation.org/openaccess/content_cvpr_2015/papers/Reininghaus_A_Stable_Multi-Scale_2015_CVPR_paper.pdf for more details.
-    :param D1: (n x 2) numpy.array encoding the (finite points of the) first diagram. Must not contain essential points (i.e. with infinite coordinate).
-    :param D2: (m x 2) numpy.array encoding the second diagram.
-    :param bandwidth: bandwidth of the Gaussian kernel with which persistence diagrams will be convolved
-    :param kernel_approx: kernel approximation class used to speed up computation. Common kernel approximations classes can be found in the scikit-learn library (such as RBFSampler for instance).
-    :returns: the persistence scale space kernel value between persistence diagrams.
-    :rtype: float 
+    
+    Parameters:
+        D1: (n x 2) numpy.array encoding the (finite points of the) first diagram. Must not contain essential points (i.e. with infinite coordinate).
+        D2: (m x 2) numpy.array encoding the second diagram.
+        bandwidth (double): bandwidth of the Gaussian kernel with which persistence diagrams will be convolved
+        kernel_approx: kernel approximation class used to speed up computation. Common kernel approximations classes can be found in the scikit-learn library (such as RBFSampler for instance).
+    
+    Returns:
+        float: the persistence scale space kernel value between persistence diagrams. 
     """
     DD1 = np.concatenate([D1, D1[:,[1,0]]], axis=0)
     DD2 = np.concatenate([D2, D2[:,[1,0]]], axis=0)
@@ -57,11 +63,14 @@ def persistence_scale_space_kernel(D1, D2, kernel_approx=None, bandwidth=1.):
 def pairwise_persistence_diagram_kernels(X, Y=None, metric="sliced_wasserstein", **kwargs):
     """
     This function computes the kernel matrix between two lists of persistence diagrams given as numpy arrays of shape (nx2).
-    :param X: first list of persistence diagrams. 
-    :param Y: second list of persistence diagrams (optional). If None, pairwise kernel values are computed from the first list only.
-    :param metric: kernel to use. It can be either a string ("sliced_wasserstein", "persistence_scale_space", "persistence_weighted_gaussian", "persistence_fisher") or a function taking two numpy arrays of shape (nx2) and (mx2) as inputs.
-    :returns: kernel matrix, i.e., numpy array of shape (num diagrams 1 x num diagrams 2)
-    :rtype: float
+
+    Parameters:    
+        X (list of n numpy arrays of shape (numx2)): first list of persistence diagrams. 
+        Y (list of m numpy arrays of shape (numx2)): second list of persistence diagrams (optional). If None, pairwise kernel values are computed from the first list only.
+        metric: kernel to use. It can be either a string ("sliced_wasserstein", "persistence_scale_space", "persistence_weighted_gaussian", "persistence_fisher") or a function taking two numpy arrays of shape (nx2) and (mx2) as inputs.
+
+    Returns: 
+        numpy array of shape (nxm): kernel matrix.
     """    
     XX = np.reshape(np.arange(len(X)), [-1,1])
     YY = None if Y is None else np.reshape(np.arange(len(Y)), [-1,1])
diff --git a/src/python/gudhi/representations/metrics.py b/src/python/gudhi/representations/metrics.py
index 4070c321..e2c30f8c 100644
--- a/src/python/gudhi/representations/metrics.py
+++ b/src/python/gudhi/representations/metrics.py
@@ -20,11 +20,14 @@ from .preprocessing import Padding
 def sliced_wasserstein_distance(D1, D2, num_directions):
     """
     This is a function for computing the sliced Wasserstein distance from two persistence diagrams. The Sliced Wasserstein distance is computed by projecting the persistence diagrams onto lines, comparing the projections with the 1-norm, and finally averaging over the lines. See http://proceedings.mlr.press/v70/carriere17a.html for more details.
-    :param D1: (n x 2) numpy.array encoding the (finite points of the) first diagram. Must not contain essential points (i.e. with infinite coordinate).
-    :param D2: (m x 2) numpy.array encoding the second diagram.
-    :param num_directions: number of lines evenly sampled from [-pi/2,pi/2] in order to approximate and speed up the distance computation.
-    :returns: the sliced Wasserstein distance between persistence diagrams.
-    :rtype: float 
+    
+    Parameters:
+        D1: (n x 2) numpy.array encoding the (finite points of the) first diagram. Must not contain essential points (i.e. with infinite coordinate).
+        D2: (m x 2) numpy.array encoding the second diagram.
+        num_directions (int): number of lines evenly sampled from [-pi/2,pi/2] in order to approximate and speed up the distance computation.
+
+    Returns: 
+        float: the sliced Wasserstein distance between persistence diagrams. 
     """
     thetas = np.linspace(-np.pi/2, np.pi/2, num=num_directions+1)[np.newaxis,:-1]
     lines = np.concatenate([np.cos(thetas), np.sin(thetas)], axis=0)
@@ -42,10 +45,13 @@ def sliced_wasserstein_distance(D1, D2, num_directions):
 def compute_persistence_diagram_projections(X, num_directions):
     """
     This is a function for projecting the points of a list of persistence diagrams (as well as their diagonal projections) onto a fixed number of lines sampled uniformly on [-pi/2, pi/2]. This function can be used as a preprocessing step in order to speed up the running time for computing all pairwise sliced Wasserstein distances / kernel values on a list of persistence diagrams. 
-    :param X: list of persistence diagrams. 
-    :param num_directions: number of lines evenly sampled from [-pi/2,pi/2] in order to approximate and speed up the distance computation.
-    :returns: list of projected persistence diagrams.
-    :rtype: float
+
+    Parameters:
+        X (list of n numpy arrays of shape (numx2)): list of persistence diagrams. 
+        num_directions (int): number of lines evenly sampled from [-pi/2,pi/2] in order to approximate and speed up the distance computation.
+
+    Returns: 
+        XX (list of n numpy arrays of shape (2*numx2)): list of projected persistence diagrams.
     """
     thetas = np.linspace(-np.pi/2, np.pi/2, num=num_directions+1)[np.newaxis,:-1]
     lines = np.concatenate([np.cos(thetas), np.sin(thetas)], axis=0)
@@ -55,10 +61,13 @@ def compute_persistence_diagram_projections(X, num_directions):
 def sliced_wasserstein_distance_on_projections(D1, D2):
     """
     This is a function for computing the sliced Wasserstein distance between two persistence diagrams that have already been projected onto some lines. It simply amounts to comparing the sorted projections with the 1-norm, and averaging over the lines. See http://proceedings.mlr.press/v70/carriere17a.html for more details.
-    :param D1: (2n x number_of_lines) numpy.array containing the n projected points of the first diagram, and the n projections of their diagonal projections.
-    :param D2: (2m x number_of_lines) numpy.array containing the m projected points of the second diagram, and the m projections of their diagonal projections.
-    :returns: the sliced Wasserstein distance between the projected persistence diagrams.
-    :rtype: float 
+
+    Parameters: 
+        D1: (2n x number_of_lines) numpy.array containing the n projected points of the first diagram, and the n projections of their diagonal projections.
+        D2: (2m x number_of_lines) numpy.array containing the m projected points of the second diagram, and the m projections of their diagonal projections.
+
+    Returns: 
+        float: the sliced Wasserstein distance between the projected persistence diagrams. 
     """
     lim1, lim2 = int(len(D1)/2), int(len(D2)/2)
     approx1, approx_diag1, approx2, approx_diag2 = D1[:lim1], D1[lim1:], D2[:lim2], D2[lim2:]
@@ -70,12 +79,15 @@ def sliced_wasserstein_distance_on_projections(D1, D2):
 def persistence_fisher_distance(D1, D2, kernel_approx=None, bandwidth=1.):
     """
     This is a function for computing the persistence Fisher distance from two 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.
-    :param D1: (n x 2) numpy.array encoding the (finite points of the) first diagram. Must not contain essential points (i.e. with infinite coordinate).
-    :param D2: (m x 2) numpy.array encoding the second diagram.
-    :param bandwidth: bandwidth of the Gaussian kernel used to turn persistence diagrams into probability distributions.
-    :param kernel_approx: kernel approximation class used to speed up computation. Common kernel approximations classes can be found in the scikit-learn library (such as RBFSampler for instance).   
-    :returns: the persistence Fisher distance between persistence diagrams.
-    :rtype: float 
+
+    Parameters: 
+        D1: (n x 2) numpy.array encoding the (finite points of the) first diagram). Must not contain essential points (i.e. with infinite coordinate).
+        D2: (m x 2) numpy.array encoding the second diagram.
+        bandwidth (float): bandwidth of the Gaussian kernel used to turn persistence diagrams into probability distributions.
+        kernel_approx: kernel approximation class used to speed up computation. Common kernel approximations classes can be found in the scikit-learn library (such as RBFSampler for instance).   
+
+    Returns: 
+        float: the persistence Fisher distance between persistence diagrams. 
     """
     projection = (1./2) * np.ones((2,2))
     diagonal_projections1 = np.matmul(D1, projection)
@@ -127,11 +139,14 @@ PAIRWISE_DISTANCE_FUNCTIONS = {
 def pairwise_persistence_diagram_distances(X, Y=None, metric="bottleneck", **kwargs):
     """
     This function computes the distance matrix between two lists of persistence diagrams given as numpy arrays of shape (nx2).
-    :param X: first list of persistence diagrams. 
-    :param Y: second list of persistence diagrams (optional). If None, pairwise distances are computed from the first list only.
-    :param metric: distance to use. It can be either a string ("sliced_wasserstein", "wasserstein", "hera_wasserstein" (Wasserstein distance computed with Hera---note that Hera is also used for the default option "wasserstein"), "pot_wasserstein" (Wasserstein distance computed with POT), "bottleneck", "persistence_fisher") or a function taking two numpy arrays of shape (nx2) and (mx2) as inputs.
-    :returns: distance matrix, i.e., numpy array of shape (num diagrams 1 x num diagrams 2)
-    :rtype: float
+
+    Parameters:
+        X (list of n numpy arrays of shape (numx2)): first list of persistence diagrams. 
+        Y (list of m numpy arrays of shape (numx2)): second list of persistence diagrams (optional). If None, pairwise distances are computed from the first list only.
+        metric: distance to use. It can be either a string ("sliced_wasserstein", "wasserstein", "hera_wasserstein" (Wasserstein distance computed with Hera---note that Hera is also used for the default option "wasserstein"), "pot_wasserstein" (Wasserstein distance computed with POT), "bottleneck", "persistence_fisher") or a function taking two numpy arrays of shape (nx2) and (mx2) as inputs.
+
+    Returns: 
+        numpy array of shape (nxm): distance matrix
     """
     XX = np.reshape(np.arange(len(X)), [-1,1])
     YY = None if Y is None else np.reshape(np.arange(len(Y)), [-1,1]) 
-- 
cgit v1.2.3


From 87311ec2d59211320e763bc9bc531858b489ff7e Mon Sep 17 00:00:00 2001
From: MathieuCarriere <mathieu.carriere3@gmail.com>
Date: Tue, 28 Apr 2020 13:28:10 -0400
Subject: added call methods + other fixes

---
 .../diagram_vectorizations_distances_kernels.py    | 98 +++++++---------------
 src/python/gudhi/representations/kernel_methods.py | 88 +++++++++++++++----
 src/python/gudhi/representations/metrics.py        | 97 +++++++++++++++++----
 src/python/gudhi/representations/preprocessing.py  | 60 +++++++++++++
 src/python/gudhi/representations/vector_methods.py | 84 +++++++++++++++++++
 5 files changed, 326 insertions(+), 101 deletions(-)

diff --git a/src/python/example/diagram_vectorizations_distances_kernels.py b/src/python/example/diagram_vectorizations_distances_kernels.py
index de22d9e7..ab7d8a16 100755
--- a/src/python/example/diagram_vectorizations_distances_kernels.py
+++ b/src/python/example/diagram_vectorizations_distances_kernels.py
@@ -9,26 +9,23 @@ from gudhi.representations import DiagramSelector, Clamping, Landscape, Silhouet
   TopologicalVector, DiagramScaler, BirthPersistenceTransform,\
   PersistenceImage, PersistenceWeightedGaussianKernel, Entropy, \
   PersistenceScaleSpaceKernel, SlicedWassersteinDistance,\
-  SlicedWassersteinKernel, BottleneckDistance, PersistenceFisherKernel
+  SlicedWassersteinKernel, BottleneckDistance, PersistenceFisherKernel, WassersteinDistance
 
-D = np.array([[0.,4.],[1.,2.],[3.,8.],[6.,8.], [0., np.inf], [5., np.inf]])
-diags = [D]
+D1 = np.array([[0.,4.],[1.,2.],[3.,8.],[6.,8.], [0., np.inf], [5., np.inf]])
 
-diags = DiagramSelector(use=True, point_type="finite").fit_transform(diags)
-diags = DiagramScaler(use=True, scalers=[([0,1], MinMaxScaler())]).fit_transform(diags)
-diags = DiagramScaler(use=True, scalers=[([1], Clamping(maximum=.9))]).fit_transform(diags)
+proc1, proc2, proc3 = DiagramSelector(use=True, point_type="finite"), DiagramScaler(use=True, scalers=[([0,1], MinMaxScaler())]), DiagramScaler(use=True, scalers=[([1], Clamping(maximum=.9))])
+D1 = proc3(proc2(proc1(D1)))
 
-D = diags[0]
-plt.scatter(D[:,0],D[:,1])
+plt.scatter(D1[:,0], D1[:,1])
 plt.plot([0.,1.],[0.,1.])
 plt.title("Test Persistence Diagram for vector methods")
 plt.show()
 
 LS = Landscape(resolution=1000)
-L = LS.fit_transform(diags)
-plt.plot(L[0][:1000])
-plt.plot(L[0][1000:2000])
-plt.plot(L[0][2000:3000])
+L = LS(D1)
+plt.plot(L[:1000])
+plt.plot(L[1000:2000])
+plt.plot(L[2000:3000])
 plt.title("Landscape")
 plt.show()
 
@@ -36,50 +33,39 @@ def pow(n):
   return lambda x: np.power(x[1]-x[0],n)
 
 SH = Silhouette(resolution=1000, weight=pow(2))
-sh = SH.fit_transform(diags)
-plt.plot(sh[0])
+plt.plot(SH(D1))
 plt.title("Silhouette")
 plt.show()
 
 BC = BettiCurve(resolution=1000)
-bc = BC.fit_transform(diags)
-plt.plot(bc[0])
+plt.plot(BC(D1))
 plt.title("Betti Curve")
 plt.show()
 
 CP = ComplexPolynomial(threshold=-1, polynomial_type="T")
-cp = CP.fit_transform(diags)
-print("Complex polynomial is " + str(cp[0,:]))
+print("Complex polynomial is " + str(CP(D1)))
 
 TV = TopologicalVector(threshold=-1)
-tv = TV.fit_transform(diags)
-print("Topological vector is " + str(tv[0,:]))
+print("Topological vector is " + str(TV(D1)))
 
 PI = PersistenceImage(bandwidth=.1, weight=lambda x: x[1], im_range=[0,1,0,1], resolution=[100,100])
-pi = PI.fit_transform(diags)
-plt.imshow(np.flip(np.reshape(pi[0], [100,100]), 0))
+plt.imshow(np.flip(np.reshape(PI(D1), [100,100]), 0))
 plt.title("Persistence Image")
 plt.show()
 
 ET = Entropy(mode="scalar")
-et = ET.fit_transform(diags)
-print("Entropy statistic is " + str(et[0,:]))
+print("Entropy statistic is " + str(ET(D1)))
 
 ET = Entropy(mode="vector", normalized=False)
-et = ET.fit_transform(diags)
-plt.plot(et[0])
+plt.plot(ET(D1))
 plt.title("Entropy function")
 plt.show()
 
-D = np.array([[1.,5.],[3.,6.],[2.,7.]])
-diags2 = [D]
+D2 = np.array([[1.,5.],[3.,6.],[2.,7.]])
+D2 = proc3(proc2(proc1(D2)))
 
-diags2 = DiagramScaler(use=True, scalers=[([0,1], MinMaxScaler())]).fit_transform(diags2)
-
-D = diags[0]
-plt.scatter(D[:,0],D[:,1])
-D = diags2[0]
-plt.scatter(D[:,0],D[:,1])
+plt.scatter(D1[:,0], D1[:,1])
+plt.scatter(D2[:,0], D2[:,1])
 plt.plot([0.,1.],[0.,1.])
 plt.title("Test Persistence Diagrams for kernel methods")
 plt.show()
@@ -88,56 +74,34 @@ def arctan(C,p):
   return lambda x: C*np.arctan(np.power(x[1], p))
 
 PWG = PersistenceWeightedGaussianKernel(bandwidth=1., kernel_approx=None, weight=arctan(1.,1.))
-X = PWG.fit(diags)
-Y = PWG.transform(diags2)
-print("PWG kernel is " + str(Y[0][0]))
+print("PWG kernel is " + str(PWG(D1, D2)))
 
 PWG = PersistenceWeightedGaussianKernel(kernel_approx=RBFSampler(gamma=1./2, n_components=100000).fit(np.ones([1,2])), weight=arctan(1.,1.))
-X = PWG.fit(diags)
-Y = PWG.transform(diags2)
-print("Approximate PWG kernel is " + str(Y[0][0]))
+print("Approximate PWG kernel is " + str(PWG(D1, D2)))
 
 PSS = PersistenceScaleSpaceKernel(bandwidth=1.)
-X = PSS.fit(diags)
-Y = PSS.transform(diags2)
-print("PSS kernel is " + str(Y[0][0]))
+print("PSS kernel is " + str(PSS(D1, D2)))
 
 PSS = PersistenceScaleSpaceKernel(kernel_approx=RBFSampler(gamma=1./2, n_components=100000).fit(np.ones([1,2])))
-X = PSS.fit(diags)
-Y = PSS.transform(diags2)
-print("Approximate PSS kernel is " + str(Y[0][0]))
+print("Approximate PSS kernel is " + str(PSS(D1, D2)))
 
 sW = SlicedWassersteinDistance(num_directions=100)
-X = sW.fit(diags)
-Y = sW.transform(diags2)
-print("SW distance is " + str(Y[0][0]))
+print("SW distance is " + str(sW(D1, D2)))
 
 SW = SlicedWassersteinKernel(num_directions=100, bandwidth=1.)
-X = SW.fit(diags)
-Y = SW.transform(diags2)
-print("SW kernel is " + str(Y[0][0]))
+print("SW kernel is " + str(SW(D1, D2)))
 
 W = WassersteinDistance(order=2, internal_p=2, mode="pot")
-X = W.fit(diags)
-Y = W.transform(diags2)
-print("Wasserstein distance (POT) is " + str(Y[0][0]))
+print("Wasserstein distance (POT) is " + str(W(D1, D2)))
 
 W = WassersteinDistance(order=2, internal_p=2, mode="hera", delta=0.0001)
-X = W.fit(diags)
-Y = W.transform(diags2)
-print("Wasserstein distance (hera) is " + str(Y[0][0]))
+print("Wasserstein distance (hera) is " + str(W(D1, D2)))
 
 W = BottleneckDistance(epsilon=.001)
-X = W.fit(diags)
-Y = W.transform(diags2)
-print("Bottleneck distance is " + str(Y[0][0]))
+print("Bottleneck distance is " + str(W(D1, D2)))
 
 PF = PersistenceFisherKernel(bandwidth_fisher=1., bandwidth=1.)
-X = PF.fit(diags)
-Y = PF.transform(diags2)
-print("PF kernel is " + str(Y[0][0]))
+print("PF kernel is " + str(PF(D1, D2)))
 
 PF = PersistenceFisherKernel(bandwidth_fisher=1., bandwidth=1., kernel_approx=RBFSampler(gamma=1./2, n_components=100000).fit(np.ones([1,2])))
-X = PF.fit(diags)
-Y = PF.transform(diags2)
-print("Approximate PF kernel is " + str(Y[0][0]))
+print("Approximate PF kernel is " + str(PF(D1, D2)))
diff --git a/src/python/gudhi/representations/kernel_methods.py b/src/python/gudhi/representations/kernel_methods.py
index 50186d63..edd1382a 100644
--- a/src/python/gudhi/representations/kernel_methods.py
+++ b/src/python/gudhi/representations/kernel_methods.py
@@ -10,14 +10,14 @@
 import numpy as np
 from sklearn.base import BaseEstimator, TransformerMixin
 from sklearn.metrics import pairwise_distances, pairwise_kernels
-from .metrics import SlicedWassersteinDistance, PersistenceFisherDistance, sklearn_wrapper, pairwise_persistence_diagram_distances, sliced_wasserstein_distance, persistence_fisher_distance
+from .metrics import SlicedWassersteinDistance, PersistenceFisherDistance, _sklearn_wrapper, pairwise_persistence_diagram_distances, _sliced_wasserstein_distance, _persistence_fisher_distance
 from .preprocessing import Padding
 
 #############################################
 # Kernel methods ############################
 #############################################
 
-def persistence_weighted_gaussian_kernel(D1, D2, weight=lambda x: 1, kernel_approx=None, bandwidth=1.):
+def _persistence_weighted_gaussian_kernel(D1, D2, weight=lambda x: 1, kernel_approx=None, bandwidth=1.):
     """
     This is a function for computing the persistence weighted Gaussian kernel value from two persistence diagrams. The persistence weighted Gaussian kernel is computed by convolving the persistence diagram points with weighted Gaussian kernels. See http://proceedings.mlr.press/v48/kusano16.html for more details.
 
@@ -25,7 +25,7 @@ def persistence_weighted_gaussian_kernel(D1, D2, weight=lambda x: 1, kernel_appr
         D1: (n x 2) numpy.array encoding the (finite points of the) first diagram. Must not contain essential points (i.e. with infinite coordinate).
         D2: (m x 2) numpy.array encoding the second diagram.
         bandwidth (double): bandwidth of the Gaussian kernel with which persistence diagrams will be convolved
-        weight: weight function for the persistence diagram points. This function must be defined on 2D points, ie lists or numpy arrays of the form [p_x,p_y].
+        weight: weight function for the persistence diagram points (default constant function, ie lambda x: 1). This function must be defined on 2D points, ie lists or numpy arrays of the form [p_x,p_y].
         kernel_approx: kernel approximation class used to speed up computation. Common kernel approximations classes can be found in the scikit-learn library (such as RBFSampler for instance).
 
     Returns:
@@ -42,7 +42,7 @@ def persistence_weighted_gaussian_kernel(D1, D2, weight=lambda x: 1, kernel_appr
         E = (1./(np.sqrt(2*np.pi)*bandwidth)) * np.exp(-np.square(pairwise_distances(D1,D2))/(2*bandwidth*bandwidth))
         return np.sum(np.multiply(W, E))
 
-def persistence_scale_space_kernel(D1, D2, kernel_approx=None, bandwidth=1.):
+def _persistence_scale_space_kernel(D1, D2, kernel_approx=None, bandwidth=1.):
     """
     This is a function for computing the persistence scale space kernel value from two persistence diagrams. The persistence scale space kernel is computed by adding the symmetric to the diagonal of each point in each persistence diagram, with negative weight, and then convolving the points with a Gaussian kernel. See https://www.cv-foundation.org/openaccess/content_cvpr_2015/papers/Reininghaus_A_Stable_Multi-Scale_2015_CVPR_paper.pdf for more details.
     
@@ -58,32 +58,32 @@ def persistence_scale_space_kernel(D1, D2, kernel_approx=None, bandwidth=1.):
     DD1 = np.concatenate([D1, D1[:,[1,0]]], axis=0)
     DD2 = np.concatenate([D2, D2[:,[1,0]]], axis=0)
     weight_pss = lambda x: 1 if x[1] >= x[0] else -1
-    return 0.5 * persistence_weighted_gaussian_kernel(DD1, DD2, weight=weight_pss, kernel_approx=kernel_approx, bandwidth=bandwidth)
+    return 0.5 * _persistence_weighted_gaussian_kernel(DD1, DD2, weight=weight_pss, kernel_approx=kernel_approx, bandwidth=bandwidth)
 
-def pairwise_persistence_diagram_kernels(X, Y=None, metric="sliced_wasserstein", **kwargs):
+def pairwise_persistence_diagram_kernels(X, Y=None, kernel="sliced_wasserstein", **kwargs):
     """
     This function computes the kernel matrix between two lists of persistence diagrams given as numpy arrays of shape (nx2).
 
     Parameters:    
         X (list of n numpy arrays of shape (numx2)): first list of persistence diagrams. 
         Y (list of m numpy arrays of shape (numx2)): second list of persistence diagrams (optional). If None, pairwise kernel values are computed from the first list only.
-        metric: kernel to use. It can be either a string ("sliced_wasserstein", "persistence_scale_space", "persistence_weighted_gaussian", "persistence_fisher") or a function taking two numpy arrays of shape (nx2) and (mx2) as inputs.
+        kernel: kernel to use. It can be either a string ("sliced_wasserstein", "persistence_scale_space", "persistence_weighted_gaussian", "persistence_fisher") or a function taking two numpy arrays of shape (nx2) and (mx2) as inputs.
 
     Returns: 
         numpy array of shape (nxm): kernel matrix.
     """    
     XX = np.reshape(np.arange(len(X)), [-1,1])
     YY = None if Y is None else np.reshape(np.arange(len(Y)), [-1,1])
-    if metric == "sliced_wasserstein":
+    if kernel == "sliced_wasserstein":
         return np.exp(-pairwise_persistence_diagram_distances(X, Y, metric="sliced_wasserstein", num_directions=kwargs["num_directions"]) / kwargs["bandwidth"])
-    elif metric == "persistence_fisher":
+    elif kernel == "persistence_fisher":
         return np.exp(-pairwise_persistence_diagram_distances(X, Y, metric="persistence_fisher", kernel_approx=kwargs["kernel_approx"], bandwidth=kwargs["bandwidth"]) / kwargs["bandwidth_fisher"])
-    elif metric == "persistence_scale_space":
-        return pairwise_kernels(XX, YY, metric=sklearn_wrapper(persistence_scale_space_kernel, X, Y, **kwargs))
-    elif metric == "persistence_weighted_gaussian":
-        return pairwise_kernels(XX, YY, metric=sklearn_wrapper(persistence_weighted_gaussian_kernel, X, Y, **kwargs))
+    elif kernel == "persistence_scale_space":
+        return pairwise_kernels(XX, YY, metric=_sklearn_wrapper(_persistence_scale_space_kernel, X, Y, **kwargs))
+    elif kernel == "persistence_weighted_gaussian":
+        return pairwise_kernels(XX, YY, metric=_sklearn_wrapper(_persistence_weighted_gaussian_kernel, X, Y, **kwargs))
     else:
-        return pairwise_kernels(XX, YY, metric=sklearn_wrapper(metric, **kwargs))
+        return pairwise_kernels(XX, YY, metric=_sklearn_wrapper(metric, **kwargs))
 
 class SlicedWassersteinKernel(BaseEstimator, TransformerMixin):
     """
@@ -121,7 +121,20 @@ class SlicedWassersteinKernel(BaseEstimator, TransformerMixin):
         Returns:
             numpy array of shape (number of diagrams in **diagrams**) x (number of diagrams in X): matrix of pairwise sliced Wasserstein kernel values.
         """
-        return pairwise_persistence_diagram_kernels(X, self.diagrams_, metric="sliced_wasserstein", bandwidth=self.bandwidth, num_directions=self.num_directions) 
+        return pairwise_persistence_diagram_kernels(X, self.diagrams_, kernel="sliced_wasserstein", bandwidth=self.bandwidth, num_directions=self.num_directions) 
+
+    def __call__(self, diag1, diag2):
+        """
+        Apply SlicedWassersteinKernel on a single pair of persistence diagrams and outputs the result.
+
+        Parameters:
+            diag1 (n x 2 numpy array): first input persistence diagram.
+            diag2 (n x 2 numpy array): second input persistence diagram.
+
+        Returns:
+            float: sliced Wasserstein kernel value.
+        """
+        return np.exp(-_sliced_wasserstein_distance(diag1, diag2, num_directions=self.num_directions)) / self.bandwidth
 
 class PersistenceWeightedGaussianKernel(BaseEstimator, TransformerMixin):
     """
@@ -160,7 +173,20 @@ class PersistenceWeightedGaussianKernel(BaseEstimator, TransformerMixin):
         Returns:
             numpy array of shape (number of diagrams in **diagrams**) x (number of diagrams in X): matrix of pairwise persistence weighted Gaussian kernel values.
         """
-        return pairwise_persistence_diagram_kernels(X, self.diagrams_, metric="persistence_weighted_gaussian", bandwidth=self.bandwidth, weight=self.weight, kernel_approx=self.kernel_approx) 
+        return pairwise_persistence_diagram_kernels(X, self.diagrams_, kernel="persistence_weighted_gaussian", bandwidth=self.bandwidth, weight=self.weight, kernel_approx=self.kernel_approx) 
+
+    def __call__(self, diag1, diag2):
+        """
+        Apply PersistenceWeightedGaussianKernel on a single pair of persistence diagrams and outputs the result.
+
+        Parameters:
+            diag1 (n x 2 numpy array): first input persistence diagram.
+            diag2 (n x 2 numpy array): second input persistence diagram.
+
+        Returns:
+            float: persistence weighted Gaussian kernel value.
+        """
+        return _persistence_weighted_gaussian_kernel(diag1, diag2, weight=self.weight, kernel_approx=self.kernel_approx, bandwidth=self.bandwidth)
 
 class PersistenceScaleSpaceKernel(BaseEstimator, TransformerMixin):
     """
@@ -197,7 +223,20 @@ class PersistenceScaleSpaceKernel(BaseEstimator, TransformerMixin):
         Returns:
             numpy array of shape (number of diagrams in **diagrams**) x (number of diagrams in X): matrix of pairwise persistence scale space kernel values.
         """
-        return pairwise_persistence_diagram_kernels(X, self.diagrams_, metric="persistence_scale_space", bandwidth=self.bandwidth, kernel_approx=self.kernel_approx)
+        return pairwise_persistence_diagram_kernels(X, self.diagrams_, kernel="persistence_scale_space", bandwidth=self.bandwidth, kernel_approx=self.kernel_approx)
+
+    def __call__(self, diag1, diag2):
+        """
+        Apply PersistenceScaleSpaceKernel on a single pair of persistence diagrams and outputs the result.
+
+        Parameters:
+            diag1 (n x 2 numpy array): first input persistence diagram.
+            diag2 (n x 2 numpy array): second input persistence diagram.
+
+        Returns:
+            float: persistence scale space kernel value.
+        """
+        return _persistence_scale_space_kernel(diag1, diag2, bandwidth=self.bandwidth, kernel_approx=self.kernel_approx)
 
 class PersistenceFisherKernel(BaseEstimator, TransformerMixin):
     """
@@ -236,5 +275,18 @@ class PersistenceFisherKernel(BaseEstimator, TransformerMixin):
         Returns:
             numpy array of shape (number of diagrams in **diagrams**) x (number of diagrams in X): matrix of pairwise persistence Fisher kernel values.
         """
-        return pairwise_persistence_diagram_kernels(X, self.diagrams_, metric="persistence_fisher", bandwidth=self.bandwidth, bandwidth_fisher=self.bandwidth_fisher, kernel_approx=self.kernel_approx)
+        return pairwise_persistence_diagram_kernels(X, self.diagrams_, kernel="persistence_fisher", bandwidth=self.bandwidth, bandwidth_fisher=self.bandwidth_fisher, kernel_approx=self.kernel_approx)
+
+    def __call__(self, diag1, diag2):
+        """
+        Apply PersistenceFisherKernel on a single pair of persistence diagrams and outputs the result.
+
+        Parameters:
+            diag1 (n x 2 numpy array): first input persistence diagram.
+            diag2 (n x 2 numpy array): second input persistence diagram.
+
+        Returns:
+            float: persistence Fisher kernel value.
+        """
+        return np.exp(-_persistence_fisher_distance(diag1, diag2, bandwidth=self.bandwidth, kernel_approx=self.kernel_approx)) / self.bandwidth_fisher
 
diff --git a/src/python/gudhi/representations/metrics.py b/src/python/gudhi/representations/metrics.py
index 59440b1a..a4bf19a6 100644
--- a/src/python/gudhi/representations/metrics.py
+++ b/src/python/gudhi/representations/metrics.py
@@ -17,7 +17,7 @@ from .preprocessing import Padding
 # Metrics ###################################
 #############################################
 
-def sliced_wasserstein_distance(D1, D2, num_directions):
+def _sliced_wasserstein_distance(D1, D2, num_directions):
     """
     This is a function for computing the sliced Wasserstein distance from two persistence diagrams. The Sliced Wasserstein distance is computed by projecting the persistence diagrams onto lines, comparing the projections with the 1-norm, and finally averaging over the lines. See http://proceedings.mlr.press/v70/carriere17a.html for more details.
     
@@ -42,7 +42,7 @@ def sliced_wasserstein_distance(D1, D2, num_directions):
     L1 = np.sum(np.abs(A-B), axis=0)
     return np.mean(L1)
 
-def compute_persistence_diagram_projections(X, num_directions):
+def _compute_persistence_diagram_projections(X, num_directions):
     """
     This is a function for projecting the points of a list of persistence diagrams (as well as their diagonal projections) onto a fixed number of lines sampled uniformly on [-pi/2, pi/2]. This function can be used as a preprocessing step in order to speed up the running time for computing all pairwise sliced Wasserstein distances / kernel values on a list of persistence diagrams. 
 
@@ -51,14 +51,14 @@ def compute_persistence_diagram_projections(X, num_directions):
         num_directions (int): number of lines evenly sampled from [-pi/2,pi/2] in order to approximate and speed up the distance computation.
 
     Returns: 
-        XX (list of n numpy arrays of shape (2*numx2)): list of projected persistence diagrams.
+        list of n numpy arrays of shape (2*numx2): list of projected persistence diagrams.
     """
     thetas = np.linspace(-np.pi/2, np.pi/2, num=num_directions+1)[np.newaxis,:-1]
     lines = np.concatenate([np.cos(thetas), np.sin(thetas)], axis=0)
     XX = [np.vstack([np.matmul(D, lines), np.matmul(np.matmul(D, .5 * np.ones((2,2))), lines)]) for D in X]
     return XX
 
-def sliced_wasserstein_distance_on_projections(D1, D2):
+def _sliced_wasserstein_distance_on_projections(D1, D2):
     """
     This is a function for computing the sliced Wasserstein distance between two persistence diagrams that have already been projected onto some lines. It simply amounts to comparing the sorted projections with the 1-norm, and averaging over the lines. See http://proceedings.mlr.press/v70/carriere17a.html for more details.
 
@@ -76,7 +76,7 @@ def sliced_wasserstein_distance_on_projections(D1, D2):
     L1 = np.sum(np.abs(A-B), axis=0)
     return np.mean(L1)
 
-def persistence_fisher_distance(D1, D2, kernel_approx=None, bandwidth=1.):
+def _persistence_fisher_distance(D1, D2, kernel_approx=None, bandwidth=1.):
     """
     This is a function for computing the persistence Fisher distance from two 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.
 
@@ -118,7 +118,7 @@ def persistence_fisher_distance(D1, D2, kernel_approx=None, bandwidth=1.):
             vectorj = vectorj/vectorj_sum
         return np.arccos(  min(np.dot(np.sqrt(vectori), np.sqrt(vectorj)), 1.)  )
 
-def sklearn_wrapper(metric, X, Y, **kwargs):
+def _sklearn_wrapper(metric, X, Y, **kwargs):
     """
     This function is a wrapper for any metric between two persistence diagrams that takes two numpy arrays of shapes (nx2) and (mx2) as arguments.
     """
@@ -133,7 +133,7 @@ def sklearn_wrapper(metric, X, Y, **kwargs):
 PAIRWISE_DISTANCE_FUNCTIONS = {
     "wasserstein": hera_wasserstein_distance,
     "hera_wasserstein": hera_wasserstein_distance,
-    "persistence_fisher": persistence_fisher_distance,
+    "persistence_fisher": _persistence_fisher_distance,
 }
 
 def pairwise_persistence_diagram_distances(X, Y=None, metric="bottleneck", **kwargs):
@@ -143,7 +143,7 @@ def pairwise_persistence_diagram_distances(X, Y=None, metric="bottleneck", **kwa
     Parameters:
         X (list of n numpy arrays of shape (numx2)): first list of persistence diagrams. 
         Y (list of m numpy arrays of shape (numx2)): second list of persistence diagrams (optional). If None, pairwise distances are computed from the first list only.
-        metric: distance to use. It can be either a string ("sliced_wasserstein", "wasserstein", "hera_wasserstein" (Wasserstein distance computed with Hera---note that Hera is also used for the default option "wasserstein"), "pot_wasserstein" (Wasserstein distance computed with POT), "bottleneck", "persistence_fisher") or a function taking two numpy arrays of shape (nx2) and (mx2) as inputs.
+        metric: distance to use. It can be either a string ("sliced_wasserstein", "wasserstein", "hera_wasserstein" (Wasserstein distance computed with Hera---note that Hera is also used for the default option "wasserstein"), "pot_wasserstein" (Wasserstein distance computed with POT), "bottleneck", "persistence_fisher") or a symmetric function taking two numpy arrays of shape (nx2) and (mx2) as inputs.
 
     Returns: 
         numpy array of shape (nxm): distance matrix
@@ -153,25 +153,25 @@ def pairwise_persistence_diagram_distances(X, Y=None, metric="bottleneck", **kwa
     if metric == "bottleneck":
         try: 
             from .. import bottleneck_distance
-            return pairwise_distances(XX, YY, metric=sklearn_wrapper(bottleneck_distance, X, Y, **kwargs))
+            return pairwise_distances(XX, YY, metric=_sklearn_wrapper(bottleneck_distance, X, Y, **kwargs))
         except ImportError:
             print("Gudhi built without CGAL")
             raise
     elif metric == "pot_wasserstein":
         try:
             from gudhi.wasserstein import wasserstein_distance as pot_wasserstein_distance
-            return pairwise_distances(XX, YY, metric=sklearn_wrapper(pot_wasserstein_distance,  X, Y, **kwargs))
+            return pairwise_distances(XX, YY, metric=_sklearn_wrapper(pot_wasserstein_distance,  X, Y, **kwargs))
         except ImportError:
             print("POT (Python Optimal Transport) is not installed. Please install POT or use metric='wasserstein' or metric='hera_wasserstein'")
             raise
     elif metric == "sliced_wasserstein":
-        Xproj = compute_persistence_diagram_projections(X, **kwargs)
-        Yproj = None if Y is None else compute_persistence_diagram_projections(Y, **kwargs)
-        return pairwise_distances(XX, YY, metric=sklearn_wrapper(sliced_wasserstein_distance_on_projections, Xproj, Yproj))
+        Xproj = _compute_persistence_diagram_projections(X, **kwargs)
+        Yproj = None if Y is None else _compute_persistence_diagram_projections(Y, **kwargs)
+        return pairwise_distances(XX, YY, metric=_sklearn_wrapper(_sliced_wasserstein_distance_on_projections, Xproj, Yproj))
     elif type(metric) == str:
-        return pairwise_distances(XX, YY, metric=sklearn_wrapper(PAIRWISE_DISTANCE_FUNCTIONS[metric], X, Y, **kwargs))
+        return pairwise_distances(XX, YY, metric=_sklearn_wrapper(PAIRWISE_DISTANCE_FUNCTIONS[metric], X, Y, **kwargs))
     else:
-        return pairwise_distances(XX, YY, metric=sklearn_wrapper(metric, X, Y, **kwargs))
+        return pairwise_distances(XX, YY, metric=_sklearn_wrapper(metric, X, Y, **kwargs))
 
 class SlicedWassersteinDistance(BaseEstimator, TransformerMixin):
     """
@@ -209,6 +209,19 @@ class SlicedWassersteinDistance(BaseEstimator, TransformerMixin):
         """
         return pairwise_persistence_diagram_distances(X, self.diagrams_, metric="sliced_wasserstein", num_directions=self.num_directions)
 
+    def __call__(self, diag1, diag2):
+        """
+        Apply SlicedWassersteinDistance on a single pair of persistence diagrams and outputs the result.
+
+        Parameters:
+            diag1 (n x 2 numpy array): first input persistence diagram.
+            diag2 (n x 2 numpy array): second input persistence diagram.
+
+        Returns:
+            float: sliced Wasserstein distance.
+        """
+        return _sliced_wasserstein_distance(diag1, diag2, num_directions=self.num_directions)
+
 class BottleneckDistance(BaseEstimator, TransformerMixin):
     """
     This is a class for computing the bottleneck distance matrix from a list of persistence diagrams. 
@@ -246,6 +259,24 @@ class BottleneckDistance(BaseEstimator, TransformerMixin):
         Xfit = pairwise_persistence_diagram_distances(X, self.diagrams_, metric="bottleneck", e=self.epsilon)
         return Xfit
 
+    def __call__(self, diag1, diag2):
+        """
+        Apply BottleneckDistance on a single pair of persistence diagrams and outputs the result.
+
+        Parameters:
+            diag1 (n x 2 numpy array): first input persistence diagram.
+            diag2 (n x 2 numpy array): second input persistence diagram.
+
+        Returns:
+            float: bottleneck distance.
+        """
+        try: 
+            from .. import bottleneck_distance
+            return bottleneck_distance(diag1, diag2, e=self.epsilon)
+        except ImportError:
+            print("Gudhi built without CGAL")
+            raise
+
 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. 
@@ -283,6 +314,19 @@ class PersistenceFisherDistance(BaseEstimator, TransformerMixin):
         """
         return pairwise_persistence_diagram_distances(X, self.diagrams_, metric="persistence_fisher", bandwidth=self.bandwidth, kernel_approx=self.kernel_approx)
 
+    def __call__(self, diag1, diag2):
+        """
+        Apply PersistenceFisherDistance on a single pair of persistence diagrams and outputs the result.
+
+        Parameters:
+            diag1 (n x 2 numpy array): first input persistence diagram.
+            diag2 (n x 2 numpy array): second input persistence diagram.
+
+        Returns:
+            float: persistence Fisher distance.
+        """
+        return _persistence_fisher_distance(diag1, diag2, bandwidth=self.bandwidth, kernel_approx=self.kernel_approx)
+
 class WassersteinDistance(BaseEstimator, TransformerMixin):
     """
     This is a class for computing the Wasserstein distance matrix from a list of persistence diagrams. 
@@ -325,5 +369,26 @@ class WassersteinDistance(BaseEstimator, TransformerMixin):
         if self.metric == "hera_wasserstein":
             Xfit = pairwise_persistence_diagram_distances(X, self.diagrams_, metric=self.metric, order=self.order, internal_p=self.internal_p, delta=self.delta)
         else:
-            Xfit = pairwise_persistence_diagram_distances(X, self.diagrams_, metric=self.metric, order=self.order, internal_p=self.internal_p)
+            Xfit = pairwise_persistence_diagram_distances(X, self.diagrams_, metric=self.metric, order=self.order, internal_p=self.internal_p, matching=False)
         return Xfit
+
+    def __call__(self, diag1, diag2):
+        """
+        Apply WassersteinDistance on a single pair of persistence diagrams and outputs the result.
+
+        Parameters:
+            diag1 (n x 2 numpy array): first input persistence diagram.
+            diag2 (n x 2 numpy array): second input persistence diagram.
+
+        Returns:
+            float: Wasserstein distance.
+        """
+        if self.metric == "hera_wasserstein":
+            return hera_wasserstein_distance(diag1, diag2, order=self.order, internal_p=self.internal_p, delta=self.delta)
+        else:
+            try:
+                from gudhi.wasserstein import wasserstein_distance as pot_wasserstein_distance
+                return pot_wasserstein_distance(diag1, diag2, order=self.order, internal_p=self.internal_p, matching=False)
+            except ImportError:
+                print("POT (Python Optimal Transport) is not installed. Please install POT or use metric='wasserstein' or metric='hera_wasserstein'")
+                raise
diff --git a/src/python/gudhi/representations/preprocessing.py b/src/python/gudhi/representations/preprocessing.py
index a39b00e4..a8545349 100644
--- a/src/python/gudhi/representations/preprocessing.py
+++ b/src/python/gudhi/representations/preprocessing.py
@@ -54,6 +54,18 @@ class BirthPersistenceTransform(BaseEstimator, TransformerMixin):
             Xfit.append(new_diag)
         return Xfit
 
+    def __call__(self, diag):
+        """
+        Apply BirthPersistenceTransform on a single persistence diagram and outputs the result.
+
+        Parameters:
+            diag (n x 2 numpy array): input persistence diagram.
+
+        Returns:
+            n x 2 numpy array: transformed persistence diagram.
+        """
+        return self.fit_transform([diag])[0]
+
 class Clamping(BaseEstimator, TransformerMixin):
     """
     This is a class for clamping values. It can be used as a parameter for the DiagramScaler class, for instance if you want to clamp abscissae or ordinates of persistence diagrams.
@@ -142,6 +154,18 @@ class DiagramScaler(BaseEstimator, TransformerMixin):
                             Xfit[i][:,I] = np.squeeze(scaler.transform(np.reshape(Xfit[i][:,I], [-1,1])))
         return Xfit
 
+    def __call__(self, diag):
+        """
+        Apply DiagramScaler on a single persistence diagram and outputs the result.
+
+        Parameters:
+            diag (n x 2 numpy array): input persistence diagram.
+
+        Returns:
+            n x 2 numpy array: transformed persistence diagram.
+        """
+        return self.fit_transform([diag])[0]
+
 class Padding(BaseEstimator, TransformerMixin):
     """
     This is a class for padding a list of persistence diagrams with dummy points, so that all persistence diagrams end up with the same number of points.
@@ -186,6 +210,18 @@ class Padding(BaseEstimator, TransformerMixin):
             Xfit = X
         return Xfit
 
+    def __call__(self, diag):
+        """
+        Apply Padding on a single persistence diagram and outputs the result.
+
+        Parameters:
+            diag (n x 2 numpy array): input persistence diagram.
+
+        Returns:
+            n x 2 numpy array: padded persistence diagram.
+        """
+        return self.fit_transform([diag])[0]
+
 class ProminentPoints(BaseEstimator, TransformerMixin):
     """
     This is a class for removing points that are close or far from the diagonal in persistence diagrams.  If persistence diagrams are n x 2 numpy arrays (i.e. persistence diagrams with ordinary features), points are ordered and thresholded by distance-to-diagonal. If persistence diagrams are n x 1 numpy arrays (i.e. persistence diagrams with essential features), points are not ordered and thresholded by first coordinate.
@@ -259,6 +295,18 @@ class ProminentPoints(BaseEstimator, TransformerMixin):
             Xfit = X
         return Xfit
 
+    def __call__(self, diag):
+        """
+        Apply ProminentPoints on a single persistence diagram and outputs the result.
+
+        Parameters:
+            diag (n x 2 numpy array): input persistence diagram.
+
+        Returns:
+            n x 2 numpy array: thresholded persistence diagram.
+        """
+        return self.fit_transform([diag])[0]
+
 class DiagramSelector(BaseEstimator, TransformerMixin):
     """
     This is a class for extracting finite or essential points in persistence diagrams.
@@ -303,3 +351,15 @@ class DiagramSelector(BaseEstimator, TransformerMixin):
         else:
             Xfit = X
         return Xfit
+
+    def __call__(self, diag):
+        """
+        Apply DiagramSelector on a single persistence diagram and outputs the result.
+
+        Parameters:
+            diag (n x 2 numpy array): input persistence diagram.
+
+        Returns:
+            n x 2 numpy array: extracted persistence diagram.
+        """
+        return self.fit_transform([diag])[0]
diff --git a/src/python/gudhi/representations/vector_methods.py b/src/python/gudhi/representations/vector_methods.py
index fe26dbe2..46fee086 100644
--- a/src/python/gudhi/representations/vector_methods.py
+++ b/src/python/gudhi/representations/vector_methods.py
@@ -81,6 +81,18 @@ class PersistenceImage(BaseEstimator, TransformerMixin):
 
         return Xfit
 
+    def __call__(self, diag):
+        """
+        Apply PersistenceImage on a single persistence diagram and outputs the result.
+
+        Parameters:
+            diag (n x 2 numpy array): input persistence diagram.
+
+        Returns:
+            numpy array with shape (number of pixels = **resolution[0]** x **resolution[1]**):: output persistence image.
+        """
+        return self.fit_transform([diag])[0,:]
+
 class Landscape(BaseEstimator, TransformerMixin):
     """
     This is a class for computing persistence landscapes from a list of persistence diagrams. A persistence landscape is a collection of 1D piecewise-linear functions computed from the rank function associated to the persistence diagram. These piecewise-linear functions are then sampled evenly on a given range and the corresponding vectors of samples are concatenated and returned. See http://jmlr.org/papers/v16/bubenik15a.html for more details.
@@ -170,6 +182,18 @@ class Landscape(BaseEstimator, TransformerMixin):
 
         return Xfit
 
+    def __call__(self, diag):
+        """
+        Apply Landscape on a single persistence diagram and outputs the result.
+
+        Parameters:
+            diag (n x 2 numpy array): input persistence diagram.
+
+        Returns:
+            numpy array with shape (number of samples = **num_landscapes** x **resolution**): output persistence landscape.
+        """
+        return self.fit_transform([diag])[0,:]
+
 class Silhouette(BaseEstimator, TransformerMixin):
     """
     This is a class for computing persistence silhouettes from a list of persistence diagrams. A persistence silhouette is computed by taking a weighted average of the collection of 1D piecewise-linear functions given by the persistence landscapes, and then by evenly sampling this average on a given range. Finally, the corresponding vector of samples is returned. See https://arxiv.org/abs/1312.0308 for more details.
@@ -248,6 +272,18 @@ class Silhouette(BaseEstimator, TransformerMixin):
 
         return Xfit 
 
+    def __call__(self, diag):
+        """
+        Apply Silhouette on a single persistence diagram and outputs the result.
+
+        Parameters:
+            diag (n x 2 numpy array): input persistence diagram.
+
+        Returns:
+            numpy array with shape (**resolution**): output persistence silhouette.
+        """
+        return self.fit_transform([diag])[0,:]
+
 class BettiCurve(BaseEstimator, TransformerMixin):
     """
     This is a class for computing Betti curves from a list of persistence diagrams. A Betti curve is a 1D piecewise-constant function obtained from the rank function. It is sampled evenly on a given range and the vector of samples is returned. See https://www.researchgate.net/publication/316604237_Time_Series_Classification_via_Topological_Data_Analysis for more details.
@@ -308,6 +344,18 @@ class BettiCurve(BaseEstimator, TransformerMixin):
 
         return Xfit
 
+    def __call__(self, diag):
+        """
+        Apply BettiCurve on a single persistence diagram and outputs the result.
+
+        Parameters:
+            diag (n x 2 numpy array): input persistence diagram.
+
+        Returns:
+            numpy array with shape (**resolution**): output Betti curve.
+        """
+        return self.fit_transform([diag])[0,:]
+
 class Entropy(BaseEstimator, TransformerMixin):
     """
     This is a class for computing persistence entropy. Persistence entropy is a statistic for persistence diagrams inspired from Shannon entropy. This statistic can also be used to compute a feature vector, called the entropy summary function. See https://arxiv.org/pdf/1803.08304.pdf for more details. Note that a previous implementation was contributed by Manuel Soriano-Trigueros.
@@ -378,6 +426,18 @@ class Entropy(BaseEstimator, TransformerMixin):
 
         return Xfit
 
+    def __call__(self, diag):
+        """
+        Apply Entropy on a single persistence diagram and outputs the result.
+
+        Parameters:
+            diag (n x 2 numpy array): input persistence diagram.
+
+        Returns:
+            numpy array with shape (1 if **mode** = "scalar" else **resolution**): output entropy.
+        """
+        return self.fit_transform([diag])[0,:]
+
 class TopologicalVector(BaseEstimator, TransformerMixin):
     """
     This is a class for computing topological vectors from a list of persistence diagrams. The topological vector associated to a persistence diagram is the sorted vector of a slight modification of the pairwise distances between the persistence diagram points. See https://diglib.eg.org/handle/10.1111/cgf12692 for more details.
@@ -431,6 +491,18 @@ class TopologicalVector(BaseEstimator, TransformerMixin):
 
         return Xfit
 
+    def __call__(self, diag):
+        """
+        Apply TopologicalVector on a single persistence diagram and outputs the result.
+
+        Parameters:
+            diag (n x 2 numpy array): input persistence diagram.
+
+        Returns:
+            numpy array with shape (**threshold**): output topological vector.
+        """
+        return self.fit_transform([diag])[0,:]
+
 class ComplexPolynomial(BaseEstimator, TransformerMixin):
     """
     This is a class for computing complex polynomials from a list of persistence diagrams. The persistence diagram points are seen as the roots of some complex polynomial, whose coefficients are returned in a complex vector. See https://link.springer.com/chapter/10.1007%2F978-3-319-23231-7_27 for more details.
@@ -490,3 +562,15 @@ class ComplexPolynomial(BaseEstimator, TransformerMixin):
             coeff = np.array(coeff[::-1])[1:]
             Xfit[d, :min(thresh, coeff.shape[0])] = coeff[:min(thresh, coeff.shape[0])]
         return Xfit
+
+    def __call__(self, diag):
+        """
+        Apply ComplexPolynomial on a single persistence diagram and outputs the result.
+
+        Parameters:
+            diag (n x 2 numpy array): input persistence diagram.
+
+        Returns:
+            numpy array with shape (**threshold**): output complex vector of coefficients.
+        """
+        return self.fit_transform([diag])[0,:]
-- 
cgit v1.2.3


From b2177e897b575e0c8d17b8ae5ed3259541a06bea Mon Sep 17 00:00:00 2001
From: MathieuCarriere <mathieu.carriere3@gmail.com>
Date: Wed, 29 Apr 2020 19:16:50 -0400
Subject: small modifs

---
 src/python/doc/representations.rst                             | 2 +-
 src/python/example/diagram_vectorizations_distances_kernels.py | 4 +++-
 src/python/gudhi/representations/kernel_methods.py             | 3 ++-
 src/python/gudhi/representations/metrics.py                    | 9 ++++-----
 4 files changed, 10 insertions(+), 8 deletions(-)

diff --git a/src/python/doc/representations.rst b/src/python/doc/representations.rst
index 11dcbcf9..041e3247 100644
--- a/src/python/doc/representations.rst
+++ b/src/python/doc/representations.rst
@@ -10,7 +10,7 @@ Representations manual
 
 This module, originally available at https://github.com/MathieuCarriere/sklearn-tda and named sklearn_tda, aims at bridging the gap between persistence diagrams and machine learning, by providing implementations of most of the vector representations for persistence diagrams in the literature, in a scikit-learn format. More specifically, it provides tools, using the scikit-learn standard interface, to compute distances and kernels on persistence diagrams, and to convert these diagrams into vectors in Euclidean space.
 
-A diagram is represented as a numpy array of shape (n,2), as can be obtained from :func:`~gudhi.SimplexTree.persistence_intervals_in_dimension` for instance. Points at infinity are represented as a numpy array of shape (n,1), storing only the birth time.
+A diagram is represented as a numpy array of shape (n,2), as can be obtained from :func:`~gudhi.SimplexTree.persistence_intervals_in_dimension` for instance. Points at infinity are represented as a numpy array of shape (n,1), storing only the birth time. The classes in this module can handle several persistence diagrams at once. In that case, the diagrams are provided as a list of numpy arrays. Note that it is not necessary for the diagrams to have the same number of points, i.e., for the corresponding arrays to have the same number of rows: all classes can handle arrays with different shapes.
 
 A small example is provided
 
diff --git a/src/python/example/diagram_vectorizations_distances_kernels.py b/src/python/example/diagram_vectorizations_distances_kernels.py
index ab7d8a16..c4a71a7a 100755
--- a/src/python/example/diagram_vectorizations_distances_kernels.py
+++ b/src/python/example/diagram_vectorizations_distances_kernels.py
@@ -13,7 +13,9 @@ from gudhi.representations import DiagramSelector, Clamping, Landscape, Silhouet
 
 D1 = np.array([[0.,4.],[1.,2.],[3.,8.],[6.,8.], [0., np.inf], [5., np.inf]])
 
-proc1, proc2, proc3 = DiagramSelector(use=True, point_type="finite"), DiagramScaler(use=True, scalers=[([0,1], MinMaxScaler())]), DiagramScaler(use=True, scalers=[([1], Clamping(maximum=.9))])
+proc1 = DiagramSelector(use=True, point_type="finite")
+proc2 = DiagramScaler(use=True, scalers=[([0,1], MinMaxScaler())])
+proc3 = DiagramScaler(use=True, scalers=[([1], Clamping(maximum=.9))])
 D1 = proc3(proc2(proc1(D1)))
 
 plt.scatter(D1[:,0], D1[:,1])
diff --git a/src/python/gudhi/representations/kernel_methods.py b/src/python/gudhi/representations/kernel_methods.py
index edd1382a..596f4f07 100644
--- a/src/python/gudhi/representations/kernel_methods.py
+++ b/src/python/gudhi/representations/kernel_methods.py
@@ -67,7 +67,8 @@ def pairwise_persistence_diagram_kernels(X, Y=None, kernel="sliced_wasserstein",
     Parameters:    
         X (list of n numpy arrays of shape (numx2)): first list of persistence diagrams. 
         Y (list of m numpy arrays of shape (numx2)): second list of persistence diagrams (optional). If None, pairwise kernel values are computed from the first list only.
-        kernel: kernel to use. It can be either a string ("sliced_wasserstein", "persistence_scale_space", "persistence_weighted_gaussian", "persistence_fisher") or a function taking two numpy arrays of shape (nx2) and (mx2) as inputs.
+        kernel: kernel to use. It can be either a string ("sliced_wasserstein", "persistence_scale_space", "persistence_weighted_gaussian", "persistence_fisher") or a function taking two numpy arrays of shape (nx2) and (mx2) as inputs. If it is a function, make sure that it is symmetric.
+        **kwargs: optional keyword parameters. Any further parameters are passed directly to the kernel function. See the docs of the various kernel classes in this module.
 
     Returns: 
         numpy array of shape (nxm): kernel matrix.
diff --git a/src/python/gudhi/representations/metrics.py b/src/python/gudhi/representations/metrics.py
index a4bf19a6..ce416fb1 100644
--- a/src/python/gudhi/representations/metrics.py
+++ b/src/python/gudhi/representations/metrics.py
@@ -32,11 +32,9 @@ def _sliced_wasserstein_distance(D1, D2, num_directions):
     thetas = np.linspace(-np.pi/2, np.pi/2, num=num_directions+1)[np.newaxis,:-1]
     lines = np.concatenate([np.cos(thetas), np.sin(thetas)], axis=0)
     approx1 = np.matmul(D1, lines)
-    diag_proj1 = (1./2) * np.ones((2,2))
-    approx_diag1 = np.matmul(np.matmul(D1, diag_proj1), lines)
+    approx_diag1 = np.matmul(np.broadcast_to(D1.sum(-1,keepdims=True)/2,(len(D1),2)), lines)
     approx2 = np.matmul(D2, lines)
-    diag_proj2 = (1./2) * np.ones((2,2))
-    approx_diag2 = np.matmul(np.matmul(D2, diag_proj2), lines)
+    approx_diag2 = np.matmul(np.broadcast_to(D2.sum(-1,keepdims=True)/2,(len(D2),2)), lines)
     A = np.sort(np.concatenate([approx1, approx_diag2], axis=0), axis=0)
     B = np.sort(np.concatenate([approx2, approx_diag1], axis=0), axis=0)
     L1 = np.sum(np.abs(A-B), axis=0)
@@ -143,7 +141,8 @@ def pairwise_persistence_diagram_distances(X, Y=None, metric="bottleneck", **kwa
     Parameters:
         X (list of n numpy arrays of shape (numx2)): first list of persistence diagrams. 
         Y (list of m numpy arrays of shape (numx2)): second list of persistence diagrams (optional). If None, pairwise distances are computed from the first list only.
-        metric: distance to use. It can be either a string ("sliced_wasserstein", "wasserstein", "hera_wasserstein" (Wasserstein distance computed with Hera---note that Hera is also used for the default option "wasserstein"), "pot_wasserstein" (Wasserstein distance computed with POT), "bottleneck", "persistence_fisher") or a symmetric function taking two numpy arrays of shape (nx2) and (mx2) as inputs.
+        metric: distance to use. It can be either a string ("sliced_wasserstein", "wasserstein", "hera_wasserstein" (Wasserstein distance computed with Hera---note that Hera is also used for the default option "wasserstein"), "pot_wasserstein" (Wasserstein distance computed with POT), "bottleneck", "persistence_fisher") or a function taking two numpy arrays of shape (nx2) and (mx2) as inputs. If it is a function, make sure that it is symmetric and that it outputs 0 if called on the same two arrays. 
+        **kwargs: optional keyword parameters. Any further parameters are passed directly to the distance function. See the docs of the various distance classes in this module.
 
     Returns: 
         numpy array of shape (nxm): distance matrix
-- 
cgit v1.2.3


From ac7917ab2cbece048e554e32cc653c14440dbcc0 Mon Sep 17 00:00:00 2001
From: Marc Glisse <marc.glisse@inria.fr>
Date: Sun, 3 May 2020 20:43:11 +0200
Subject: Fewer copies and no GIL for hera

Now the input arrays are not copied as long as they use a float64 data
type, even if they are not contiguous. That's not important here, but I
wanted an example of how to do it.
More importantly, no need to hold the GIL.
I was too lazy to benchmark to see if that changed anything...
---
 src/python/gudhi/hera.cc | 28 ++++++++++++++++++++--------
 1 file changed, 20 insertions(+), 8 deletions(-)

diff --git a/src/python/gudhi/hera.cc b/src/python/gudhi/hera.cc
index 0d562b4c..50d49c77 100644
--- a/src/python/gudhi/hera.cc
+++ b/src/python/gudhi/hera.cc
@@ -11,14 +11,24 @@
 #include <pybind11/pybind11.h>
 #include <pybind11/numpy.h>
 
-#include <boost/range/iterator_range.hpp>
+#include <boost/range/counting_range.hpp>
+#include <boost/range/adaptor/transformed.hpp>
 
 #include <wasserstein.h> // Hera
 
-#include <array>
+#include <utility>
 
 namespace py = pybind11;
-typedef py::array_t<double, py::array::c_style | py::array::forcecast> Dgm;
+typedef py::array_t<double> Dgm;
+
+// Get m[i,0] and m[i,1] as a pair
+auto pairify(void* p, ssize_t h, ssize_t w) {
+  return [=](ssize_t i){
+    char* birth = (char*)p + i * h;
+    char* death = birth + w;
+    return std::make_pair(*(double*)birth, *(double*)death);
+  };
+}
 
 double wasserstein_distance(
     Dgm d1, Dgm d2,
@@ -27,16 +37,18 @@ double wasserstein_distance(
 {
   py::buffer_info buf1 = d1.request();
   py::buffer_info buf2 = d2.request();
+
+  py::gil_scoped_release release;
+
   // shape (n,2) or (0) for empty
   if((buf1.ndim!=2 || buf1.shape[1]!=2) && (buf1.ndim!=1 || buf1.shape[0]!=0))
     throw std::runtime_error("Diagram 1 must be an array of size n x 2");
   if((buf2.ndim!=2 || buf2.shape[1]!=2) && (buf2.ndim!=1 || buf2.shape[0]!=0))
     throw std::runtime_error("Diagram 2 must be an array of size n x 2");
-  typedef std::array<double, 2> Point;
-  auto p1 = (Point*)buf1.ptr;
-  auto p2 = (Point*)buf2.ptr;
-  auto diag1 = boost::make_iterator_range(p1, p1+buf1.shape[0]);
-  auto diag2 = boost::make_iterator_range(p2, p2+buf2.shape[0]);
+  auto cnt1 = boost::counting_range<ssize_t>(0, buf1.shape[0]);
+  auto diag1 = boost::adaptors::transform(cnt1, pairify(buf1.ptr, buf1.strides[0], buf1.strides[1]));
+  auto cnt2 = boost::counting_range<ssize_t>(0, buf2.shape[0]);
+  auto diag2 = boost::adaptors::transform(cnt2, pairify(buf2.ptr, buf2.strides[0], buf2.strides[1]));
 
   hera::AuctionParams<double> params;
   params.wasserstein_power = wasserstein_power;
-- 
cgit v1.2.3


From d2a9aed9ada419b7715a77322ad17ddf3535d133 Mon Sep 17 00:00:00 2001
From: ROUVREAU Vincent <vincent.rouvreau@inria.fr>
Date: Mon, 4 May 2020 19:23:40 +0200
Subject: Try to build with conda as brew fails

---
 azure-pipelines.yml | 35 ++++++++++++++++++-----------------
 1 file changed, 18 insertions(+), 17 deletions(-)

diff --git a/azure-pipelines.yml b/azure-pipelines.yml
index 95b15db2..fccb7d57 100644
--- a/azure-pipelines.yml
+++ b/azure-pipelines.yml
@@ -4,35 +4,36 @@ jobs:
     displayName: "Build and test"
     timeoutInMinutes: 0
     cancelTimeoutInMinutes: 60
-  
+    pool:
+      vmImage: macOS-10.14
     strategy:
       matrix:
-        macOSrelease:
-          imageName: 'macos-10.14'
-          CMakeBuildType: Release
-          customInstallation: 'brew update && brew install graphviz doxygen boost eigen gmp mpfr tbb cgal'
+        Python36:
+          python.version: '3.6'
+        Python37:
+          python.version: '3.7'
+        Python38:
+          python.version: '3.8'
 
-    pool:
-      vmImage: $(imageName)
-  
     steps:
-    - task: UsePythonVersion@0
-      inputs:
-        versionSpec: '3.7'
-        architecture: 'x64'
+    - bash: echo "##vso[task.prependpath]$CONDA/bin"
+      displayName: Add conda to PATH
+  
+    - bash: conda create --yes --quiet --name gudhi_build_env
+      displayName: Create Anaconda environment
   
-    - script: |
-        $(customInstallation)
+    - bash: |
+        source activate gudhi_build_env
+        conda install --yes --quiet --name gudhi_build_env python=$PYTHON_VERSION
         git submodule update --init
-        python -m pip install --upgrade pip
         python -m pip install --user -r .github/build-requirements.txt
         python -m pip install --user -r .github/test-requirements.txt
       displayName: 'Install build dependencies'
-    - script: |
+    - bash: |
         mkdir build
         cd build
         cmake -DCMAKE_BUILD_TYPE:STRING=$(CMakeBuildType) -DWITH_GUDHI_TEST=ON -DWITH_GUDHI_UTILITIES=ON -DWITH_GUDHI_PYTHON=ON -DPython_ADDITIONAL_VERSIONS=3 ..
         make
         make doxygen
-        ctest -j 8 --output-on-failure -E sphinx # remove sphinx build as it fails
+        ctest -j 8 --output-on-failure # -E sphinx remove sphinx build as it fails
       displayName: 'Build, test and documentation generation'
-- 
cgit v1.2.3


From 5d03351f22f2511e3f5159f19f54b21bf2a04d61 Mon Sep 17 00:00:00 2001
From: ROUVREAU Vincent <vincent.rouvreau@inria.fr>
Date: Mon, 4 May 2020 19:29:02 +0200
Subject: sudo ?

---
 azure-pipelines.yml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/azure-pipelines.yml b/azure-pipelines.yml
index fccb7d57..b50bd91a 100644
--- a/azure-pipelines.yml
+++ b/azure-pipelines.yml
@@ -19,7 +19,7 @@ jobs:
     - bash: echo "##vso[task.prependpath]$CONDA/bin"
       displayName: Add conda to PATH
   
-    - bash: conda create --yes --quiet --name gudhi_build_env
+    - bash: sudo conda create --yes --quiet --name gudhi_build_env
       displayName: Create Anaconda environment
   
     - bash: |
-- 
cgit v1.2.3


From 7bd1941307033193da4c1cfcef873e69ca7f68f3 Mon Sep 17 00:00:00 2001
From: ROUVREAU Vincent <vincent.rouvreau@inria.fr>
Date: Mon, 4 May 2020 19:30:55 +0200
Subject: sudo ?

---
 azure-pipelines.yml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/azure-pipelines.yml b/azure-pipelines.yml
index b50bd91a..0fea11f6 100644
--- a/azure-pipelines.yml
+++ b/azure-pipelines.yml
@@ -24,7 +24,7 @@ jobs:
   
     - bash: |
         source activate gudhi_build_env
-        conda install --yes --quiet --name gudhi_build_env python=$PYTHON_VERSION
+        sudo conda install --yes --quiet --name gudhi_build_env python=$PYTHON_VERSION
         git submodule update --init
         python -m pip install --user -r .github/build-requirements.txt
         python -m pip install --user -r .github/test-requirements.txt
-- 
cgit v1.2.3


From 62139c92181b7f405ce0e36ef6b46777cee85b34 Mon Sep 17 00:00:00 2001
From: ROUVREAU Vincent <vincent.rouvreau@inria.fr>
Date: Mon, 4 May 2020 22:22:26 +0200
Subject: Add conda build requirements

---
 azure-pipelines.yml | 1 +
 1 file changed, 1 insertion(+)

diff --git a/azure-pipelines.yml b/azure-pipelines.yml
index 0fea11f6..97c84136 100644
--- a/azure-pipelines.yml
+++ b/azure-pipelines.yml
@@ -25,6 +25,7 @@ jobs:
     - bash: |
         source activate gudhi_build_env
         sudo conda install --yes --quiet --name gudhi_build_env python=$PYTHON_VERSION
+        sudo conda install --yes -c conda-forge doxygen eigen boost-cpp=1.70.0 cgal-cpp>=5.0
         git submodule update --init
         python -m pip install --user -r .github/build-requirements.txt
         python -m pip install --user -r .github/test-requirements.txt
-- 
cgit v1.2.3


From b880228fb423aeb3d662416fbb477d3ced100e08 Mon Sep 17 00:00:00 2001
From: ROUVREAU Vincent <vincent.rouvreau@inria.fr>
Date: Mon, 4 May 2020 22:31:13 +0200
Subject: Need to activate conda env to build

---
 azure-pipelines.yml | 11 ++++++-----
 1 file changed, 6 insertions(+), 5 deletions(-)

diff --git a/azure-pipelines.yml b/azure-pipelines.yml
index 97c84136..2fcff411 100644
--- a/azure-pipelines.yml
+++ b/azure-pipelines.yml
@@ -10,10 +10,10 @@ jobs:
       matrix:
         Python36:
           python.version: '3.6'
-        Python37:
-          python.version: '3.7'
-        Python38:
-          python.version: '3.8'
+        #Python37:
+        #  python.version: '3.7'
+        #Python38:
+        #  python.version: '3.8'
 
     steps:
     - bash: echo "##vso[task.prependpath]$CONDA/bin"
@@ -26,11 +26,12 @@ jobs:
         source activate gudhi_build_env
         sudo conda install --yes --quiet --name gudhi_build_env python=$PYTHON_VERSION
         sudo conda install --yes -c conda-forge doxygen eigen boost-cpp=1.70.0 cgal-cpp>=5.0
-        git submodule update --init
         python -m pip install --user -r .github/build-requirements.txt
         python -m pip install --user -r .github/test-requirements.txt
       displayName: 'Install build dependencies'
     - bash: |
+        source activate gudhi_build_env
+        git submodule update --init
         mkdir build
         cd build
         cmake -DCMAKE_BUILD_TYPE:STRING=$(CMakeBuildType) -DWITH_GUDHI_TEST=ON -DWITH_GUDHI_UTILITIES=ON -DWITH_GUDHI_PYTHON=ON -DPython_ADDITIONAL_VERSIONS=3 ..
-- 
cgit v1.2.3


From 71d958891cc638b26541ca5cf6c569b43332d2b6 Mon Sep 17 00:00:00 2001
From: ROUVREAU Vincent <vincent.rouvreau@inria.fr>
Date: Mon, 4 May 2020 22:39:10 +0200
Subject: conda update and release cmake version

---
 azure-pipelines.yml | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/azure-pipelines.yml b/azure-pipelines.yml
index 2fcff411..b3b0ea7f 100644
--- a/azure-pipelines.yml
+++ b/azure-pipelines.yml
@@ -10,6 +10,7 @@ jobs:
       matrix:
         Python36:
           python.version: '3.6'
+          CMakeBuildType: Release
         #Python37:
         #  python.version: '3.7'
         #Python38:
@@ -25,6 +26,7 @@ jobs:
     - bash: |
         source activate gudhi_build_env
         sudo conda install --yes --quiet --name gudhi_build_env python=$PYTHON_VERSION
+        sudo conda update --yes --quiet -n base -c defaults conda
         sudo conda install --yes -c conda-forge doxygen eigen boost-cpp=1.70.0 cgal-cpp>=5.0
         python -m pip install --user -r .github/build-requirements.txt
         python -m pip install --user -r .github/test-requirements.txt
-- 
cgit v1.2.3


From 03b8322e9ded09cc879867008d32baa3a91a45e5 Mon Sep 17 00:00:00 2001
From: ROUVREAU Vincent <vincent.rouvreau@inria.fr>
Date: Tue, 5 May 2020 07:05:45 +0200
Subject: brew install cgal & Cie instead of conda install because of link
 issue

---
 azure-pipelines.yml | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/azure-pipelines.yml b/azure-pipelines.yml
index b3b0ea7f..3ab2f112 100644
--- a/azure-pipelines.yml
+++ b/azure-pipelines.yml
@@ -11,6 +11,7 @@ jobs:
         Python36:
           python.version: '3.6'
           CMakeBuildType: Release
+          customInstallation: 'brew update && brew install graphviz doxygen boost eigen gmp mpfr tbb cgal'
         #Python37:
         #  python.version: '3.7'
         #Python38:
@@ -26,10 +27,9 @@ jobs:
     - bash: |
         source activate gudhi_build_env
         sudo conda install --yes --quiet --name gudhi_build_env python=$PYTHON_VERSION
-        sudo conda update --yes --quiet -n base -c defaults conda
-        sudo conda install --yes -c conda-forge doxygen eigen boost-cpp=1.70.0 cgal-cpp>=5.0
         python -m pip install --user -r .github/build-requirements.txt
         python -m pip install --user -r .github/test-requirements.txt
+        $(customInstallation)
       displayName: 'Install build dependencies'
     - bash: |
         source activate gudhi_build_env
-- 
cgit v1.2.3


From 8da9158e9a2ffb128eb1b5b05d4e8574ff70d771 Mon Sep 17 00:00:00 2001
From: ROUVREAU Vincent <vincent.rouvreau@inria.fr>
Date: Tue, 5 May 2020 07:48:16 +0200
Subject: Remove sphinx test and matrix

---
 azure-pipelines.yml | 22 ++++++++--------------
 1 file changed, 8 insertions(+), 14 deletions(-)

diff --git a/azure-pipelines.yml b/azure-pipelines.yml
index 3ab2f112..7b5334a7 100644
--- a/azure-pipelines.yml
+++ b/azure-pipelines.yml
@@ -6,16 +6,10 @@ jobs:
     cancelTimeoutInMinutes: 60
     pool:
       vmImage: macOS-10.14
-    strategy:
-      matrix:
-        Python36:
-          python.version: '3.6'
-          CMakeBuildType: Release
-          customInstallation: 'brew update && brew install graphviz doxygen boost eigen gmp mpfr tbb cgal'
-        #Python37:
-        #  python.version: '3.7'
-        #Python38:
-        #  python.version: '3.8'
+    variables:
+      pythonVersion: '3.6'
+      cmakeBuildType: Release
+      customInstallation: 'brew update && brew install graphviz doxygen boost eigen gmp mpfr tbb cgal'
 
     steps:
     - bash: echo "##vso[task.prependpath]$CONDA/bin"
@@ -26,7 +20,7 @@ jobs:
   
     - bash: |
         source activate gudhi_build_env
-        sudo conda install --yes --quiet --name gudhi_build_env python=$PYTHON_VERSION
+        sudo conda install --yes --quiet --name gudhi_build_env python=$(pythonVersion)
         python -m pip install --user -r .github/build-requirements.txt
         python -m pip install --user -r .github/test-requirements.txt
         $(customInstallation)
@@ -36,8 +30,8 @@ jobs:
         git submodule update --init
         mkdir build
         cd build
-        cmake -DCMAKE_BUILD_TYPE:STRING=$(CMakeBuildType) -DWITH_GUDHI_TEST=ON -DWITH_GUDHI_UTILITIES=ON -DWITH_GUDHI_PYTHON=ON -DPython_ADDITIONAL_VERSIONS=3 ..
-        make
+        cmake -DCMAKE_BUILD_TYPE:STRING=$(cmakeBuildType) -DWITH_GUDHI_TEST=ON -DWITH_GUDHI_UTILITIES=ON -DWITH_GUDHI_PYTHON=ON -DPython_ADDITIONAL_VERSIONS=3 ..
+        make -j 4
         make doxygen
-        ctest -j 8 --output-on-failure # -E sphinx remove sphinx build as it fails
+        ctest -j 4 --output-on-failure -E sphinx # remove sphinx build as it fails
       displayName: 'Build, test and documentation generation'
-- 
cgit v1.2.3


From 99549c20e9173b536ac816ab683bc13025f182a2 Mon Sep 17 00:00:00 2001
From: Marc Glisse <marc.glisse@inria.fr>
Date: Tue, 5 May 2020 11:07:53 +0200
Subject: fix use of threads and n_jobs in Parallel

---
 src/python/gudhi/point_cloud/knn.py | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/src/python/gudhi/point_cloud/knn.py b/src/python/gudhi/point_cloud/knn.py
index 07553d6d..34e80b5d 100644
--- a/src/python/gudhi/point_cloud/knn.py
+++ b/src/python/gudhi/point_cloud/knn.py
@@ -200,8 +200,8 @@ class KNearestNeighbors:
                     from joblib import Parallel, delayed, effective_n_jobs
                     from sklearn.utils import gen_even_slices
 
-                    slices = gen_even_slices(len(X), effective_n_jobs(-1))
-                    parallel = Parallel(backend="threading", n_jobs=-1)
+                    slices = gen_even_slices(len(X), effective_n_jobs(n_jobs))
+                    parallel = Parallel(prefer="threads", n_jobs=n_jobs)
                     if self.params.get("sort_results", True):
 
                         def func(M):
@@ -242,8 +242,8 @@ class KNearestNeighbors:
 
                     else:
                         func = lambda M: numpy.partition(M, k - 1)[:, 0:k]
-                    slices = gen_even_slices(len(X), effective_n_jobs(-1))
-                    parallel = Parallel(backend="threading", n_jobs=-1)
+                    slices = gen_even_slices(len(X), effective_n_jobs(n_jobs))
+                    parallel = Parallel(prefer="threads", n_jobs=n_jobs)
                     distances = numpy.concatenate(parallel(delayed(func)(X[s]) for s in slices))
                 return distances
             return None
-- 
cgit v1.2.3


From dac92c5ae9da6aa21fdcd261737e08d6898dbbdc Mon Sep 17 00:00:00 2001
From: Marc Glisse <marc.glisse@inria.fr>
Date: Wed, 6 May 2020 12:54:21 +0200
Subject: Avoid reading outside of allocated region

The result was unused, but better be safe.
---
 src/python/gudhi/hera.cc | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/src/python/gudhi/hera.cc b/src/python/gudhi/hera.cc
index 50d49c77..63bbb075 100644
--- a/src/python/gudhi/hera.cc
+++ b/src/python/gudhi/hera.cc
@@ -45,10 +45,12 @@ double wasserstein_distance(
     throw std::runtime_error("Diagram 1 must be an array of size n x 2");
   if((buf2.ndim!=2 || buf2.shape[1]!=2) && (buf2.ndim!=1 || buf2.shape[0]!=0))
     throw std::runtime_error("Diagram 2 must be an array of size n x 2");
+  ssize_t stride11 = buf1.ndim == 2 ? buf1.strides[1] : 0;
+  ssize_t stride21 = buf2.ndim == 2 ? buf2.strides[1] : 0;
   auto cnt1 = boost::counting_range<ssize_t>(0, buf1.shape[0]);
-  auto diag1 = boost::adaptors::transform(cnt1, pairify(buf1.ptr, buf1.strides[0], buf1.strides[1]));
+  auto diag1 = boost::adaptors::transform(cnt1, pairify(buf1.ptr, buf1.strides[0], stride11));
   auto cnt2 = boost::counting_range<ssize_t>(0, buf2.shape[0]);
-  auto diag2 = boost::adaptors::transform(cnt2, pairify(buf2.ptr, buf2.strides[0], buf2.strides[1]));
+  auto diag2 = boost::adaptors::transform(cnt2, pairify(buf2.ptr, buf2.strides[0], stride21));
 
   hera::AuctionParams<double> params;
   params.wasserstein_power = wasserstein_power;
-- 
cgit v1.2.3


From 5c5e2c3075235079fda94fc6a159cc5275f85a0c Mon Sep 17 00:00:00 2001
From: Marc Glisse <marc.glisse@inria.fr>
Date: Wed, 6 May 2020 14:13:14 +0200
Subject: Refactor the numpy -> C++ range conversion

If we want to reuse it for bottleneck...
---
 src/python/gudhi/hera.cc | 31 ++++++++++++++++---------------
 1 file changed, 16 insertions(+), 15 deletions(-)

diff --git a/src/python/gudhi/hera.cc b/src/python/gudhi/hera.cc
index 63bbb075..5aec1806 100644
--- a/src/python/gudhi/hera.cc
+++ b/src/python/gudhi/hera.cc
@@ -22,7 +22,7 @@ namespace py = pybind11;
 typedef py::array_t<double> Dgm;
 
 // Get m[i,0] and m[i,1] as a pair
-auto pairify(void* p, ssize_t h, ssize_t w) {
+static auto pairify(void* p, ssize_t h, ssize_t w) {
   return [=](ssize_t i){
     char* birth = (char*)p + i * h;
     char* death = birth + w;
@@ -30,28 +30,29 @@ auto pairify(void* p, ssize_t h, ssize_t w) {
   };
 }
 
+inline auto numpy_to_range_of_pairs(py::array_t<double> dgm) {
+  py::buffer_info buf = dgm.request();
+  // shape (n,2) or (0) for empty
+  if((buf.ndim!=2 || buf.shape[1]!=2) && (buf.ndim!=1 || buf.shape[0]!=0))
+    throw std::runtime_error("Diagram must be an array of size n x 2");
+  // In the case of shape (0), avoid reading non-existing strides[1] even if we won't use it.
+  ssize_t stride1 = buf.ndim == 2 ? buf.strides[1] : 0;
+  auto cnt = boost::counting_range<ssize_t>(0, buf.shape[0]);
+  return boost::adaptors::transform(cnt, pairify(buf.ptr, buf.strides[0], stride1));
+  // Be careful that the returned range cannot contain references to dead temporaries.
+}
+
 double wasserstein_distance(
     Dgm d1, Dgm d2,
     double wasserstein_power, double internal_p,
     double delta)
 {
-  py::buffer_info buf1 = d1.request();
-  py::buffer_info buf2 = d2.request();
+  // I *think* the call to request() has to be before releasing the GIL.
+  auto diag1 = numpy_to_range_of_pairs(d1);
+  auto diag2 = numpy_to_range_of_pairs(d2);
 
   py::gil_scoped_release release;
 
-  // shape (n,2) or (0) for empty
-  if((buf1.ndim!=2 || buf1.shape[1]!=2) && (buf1.ndim!=1 || buf1.shape[0]!=0))
-    throw std::runtime_error("Diagram 1 must be an array of size n x 2");
-  if((buf2.ndim!=2 || buf2.shape[1]!=2) && (buf2.ndim!=1 || buf2.shape[0]!=0))
-    throw std::runtime_error("Diagram 2 must be an array of size n x 2");
-  ssize_t stride11 = buf1.ndim == 2 ? buf1.strides[1] : 0;
-  ssize_t stride21 = buf2.ndim == 2 ? buf2.strides[1] : 0;
-  auto cnt1 = boost::counting_range<ssize_t>(0, buf1.shape[0]);
-  auto diag1 = boost::adaptors::transform(cnt1, pairify(buf1.ptr, buf1.strides[0], stride11));
-  auto cnt2 = boost::counting_range<ssize_t>(0, buf2.shape[0]);
-  auto diag2 = boost::adaptors::transform(cnt2, pairify(buf2.ptr, buf2.strides[0], stride21));
-
   hera::AuctionParams<double> params;
   params.wasserstein_power = wasserstein_power;
   // hera encodes infinity as -1...
-- 
cgit v1.2.3


From 47e5ac79af3a354358515c0213b28848f878fde6 Mon Sep 17 00:00:00 2001
From: Marc Glisse <marc.glisse@inria.fr>
Date: Wed, 6 May 2020 22:59:36 +0200
Subject: Reimplement the bottleneck python wrapper with pybind11

---
 src/python/CMakeLists.txt                   | 33 ++++++++++---------
 src/python/gudhi/bottleneck.cc              | 51 +++++++++++++++++++++++++++++
 src/python/gudhi/bottleneck.pyx             | 48 ---------------------------
 src/python/gudhi/hera.cc                    | 32 +-----------------
 src/python/include/pybind11_diagram_utils.h | 39 ++++++++++++++++++++++
 src/python/setup.py.in                      | 19 +++++++++--
 6 files changed, 125 insertions(+), 97 deletions(-)
 create mode 100644 src/python/gudhi/bottleneck.cc
 delete mode 100644 src/python/gudhi/bottleneck.pyx
 create mode 100644 src/python/include/pybind11_diagram_utils.h

diff --git a/src/python/CMakeLists.txt b/src/python/CMakeLists.txt
index d712e189..976a8b52 100644
--- a/src/python/CMakeLists.txt
+++ b/src/python/CMakeLists.txt
@@ -34,6 +34,7 @@ endfunction( add_gudhi_debug_info )
 if(PYTHONINTERP_FOUND)
   if(PYBIND11_FOUND)
     add_gudhi_debug_info("Pybind11 version ${PYBIND11_VERSION}")
+    set(GUDHI_PYTHON_MODULES "${GUDHI_PYTHON_MODULES}'bottleneck', ")
     set(GUDHI_PYTHON_MODULES_EXTRA "${GUDHI_PYTHON_MODULES_EXTRA}'hera', ")
   endif()
   if(CYTHON_FOUND)
@@ -46,7 +47,6 @@ if(PYTHONINTERP_FOUND)
     set(GUDHI_PYTHON_MODULES "${GUDHI_PYTHON_MODULES}'reader_utils', ")
     set(GUDHI_PYTHON_MODULES "${GUDHI_PYTHON_MODULES}'witness_complex', ")
     set(GUDHI_PYTHON_MODULES "${GUDHI_PYTHON_MODULES}'strong_witness_complex', ")
-    set(GUDHI_PYTHON_MODULES "${GUDHI_PYTHON_MODULES}'bottleneck', ")
     set(GUDHI_PYTHON_MODULES "${GUDHI_PYTHON_MODULES}'nerve_gic', ")
     set(GUDHI_PYTHON_MODULES "${GUDHI_PYTHON_MODULES}'subsampling', ")
     set(GUDHI_PYTHON_MODULES "${GUDHI_PYTHON_MODULES}'tangential_complex', ")
@@ -120,24 +120,25 @@ if(PYTHONINTERP_FOUND)
       set(GUDHI_PYTHON_EXTRA_COMPILE_ARGS "${GUDHI_PYTHON_EXTRA_COMPILE_ARGS}'-DCGAL_EIGEN3_ENABLED', ")
     endif (EIGEN3_FOUND)
 
-    set(GUDHI_PYTHON_MODULES_TO_COMPILE "${GUDHI_PYTHON_MODULES_TO_COMPILE}'off_reader', ")
-    set(GUDHI_PYTHON_MODULES_TO_COMPILE "${GUDHI_PYTHON_MODULES_TO_COMPILE}'simplex_tree', ")
-    set(GUDHI_PYTHON_MODULES_TO_COMPILE "${GUDHI_PYTHON_MODULES_TO_COMPILE}'rips_complex', ")
-    set(GUDHI_PYTHON_MODULES_TO_COMPILE "${GUDHI_PYTHON_MODULES_TO_COMPILE}'cubical_complex', ")
-    set(GUDHI_PYTHON_MODULES_TO_COMPILE "${GUDHI_PYTHON_MODULES_TO_COMPILE}'periodic_cubical_complex', ")
-    set(GUDHI_PYTHON_MODULES_TO_COMPILE "${GUDHI_PYTHON_MODULES_TO_COMPILE}'reader_utils', ")
-    set(GUDHI_PYTHON_MODULES_TO_COMPILE "${GUDHI_PYTHON_MODULES_TO_COMPILE}'witness_complex', ")
-    set(GUDHI_PYTHON_MODULES_TO_COMPILE "${GUDHI_PYTHON_MODULES_TO_COMPILE}'strong_witness_complex', ")
+    set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'off_reader', ")
+    set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'simplex_tree', ")
+    set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'rips_complex', ")
+    set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'cubical_complex', ")
+    set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'periodic_cubical_complex', ")
+    set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'reader_utils', ")
+    set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'witness_complex', ")
+    set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'strong_witness_complex', ")
+    set(GUDHI_PYBIND11_MODULES "${GUDHI_PYBIND11_MODULES}'hera', ")
     if (NOT CGAL_VERSION VERSION_LESS 4.11.0)
-      set(GUDHI_PYTHON_MODULES_TO_COMPILE "${GUDHI_PYTHON_MODULES_TO_COMPILE}'bottleneck', ")
-      set(GUDHI_PYTHON_MODULES_TO_COMPILE "${GUDHI_PYTHON_MODULES_TO_COMPILE}'nerve_gic', ")
+      set(GUDHI_PYBIND11_MODULES "${GUDHI_PYBIND11_MODULES}'bottleneck', ")
+      set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'nerve_gic', ")
     endif ()
     if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0)
-      set(GUDHI_PYTHON_MODULES_TO_COMPILE "${GUDHI_PYTHON_MODULES_TO_COMPILE}'alpha_complex', ")
-      set(GUDHI_PYTHON_MODULES_TO_COMPILE "${GUDHI_PYTHON_MODULES_TO_COMPILE}'subsampling', ")
-      set(GUDHI_PYTHON_MODULES_TO_COMPILE "${GUDHI_PYTHON_MODULES_TO_COMPILE}'tangential_complex', ")
-      set(GUDHI_PYTHON_MODULES_TO_COMPILE "${GUDHI_PYTHON_MODULES_TO_COMPILE}'euclidean_witness_complex', ")
-      set(GUDHI_PYTHON_MODULES_TO_COMPILE "${GUDHI_PYTHON_MODULES_TO_COMPILE}'euclidean_strong_witness_complex', ")
+      set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'alpha_complex', ")
+      set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'subsampling', ")
+      set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'tangential_complex', ")
+      set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'euclidean_witness_complex', ")
+      set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'euclidean_strong_witness_complex', ")
     endif ()
 
     if(CGAL_FOUND)
diff --git a/src/python/gudhi/bottleneck.cc b/src/python/gudhi/bottleneck.cc
new file mode 100644
index 00000000..577e5e0b
--- /dev/null
+++ b/src/python/gudhi/bottleneck.cc
@@ -0,0 +1,51 @@
+/*    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):       Marc Glisse
+ *
+ *    Copyright (C) 2020 Inria
+ *
+ *    Modification(s):
+ *      - YYYY/MM Author: Description of the modification
+ */
+
+#include <gudhi/Bottleneck.h>
+
+#include <pybind11_diagram_utils.h>
+
+double bottleneck(Dgm d1, Dgm d2, double epsilon)
+{
+  // I *think* the call to request() has to be before releasing the GIL.
+  auto diag1 = numpy_to_range_of_pairs(d1);
+  auto diag2 = numpy_to_range_of_pairs(d2);
+
+  py::gil_scoped_release release;
+
+  return Gudhi::persistence_diagram::bottleneck_distance(diag1, diag2, epsilon);
+}
+
+PYBIND11_MODULE(bottleneck, m) {
+      m.attr("__license__") = "GPL v3";
+      m.def("bottleneck_distance", &bottleneck,
+          py::arg("diagram_1"), py::arg("diagram_2"),
+          py::arg("e") = (std::numeric_limits<double>::min)(),
+          R"pbdoc(
+    This function returns the point corresponding to a given vertex.
+
+    :param diagram_1: The first diagram.
+    :type diagram_1: vector[pair[double, double]]
+    :param diagram_2: The second diagram.
+    :type diagram_2: vector[pair[double, double]]
+    :param e: If `e` is 0, this uses an expensive algorithm to compute the
+        exact distance.
+        If `e` is not 0, it asks for an additive `e`-approximation, and
+        currently also allows a small multiplicative error (the last 2 or 3
+        bits of the mantissa may be wrong). This version of the algorithm takes
+        advantage of the limited precision of `double` and is usually a lot
+        faster to compute, whatever the value of `e`.
+
+        Thus, by default, `e` is the smallest positive double.
+    :type e: float
+    :rtype: float
+    :returns: the bottleneck distance.
+    )pbdoc");
+}
diff --git a/src/python/gudhi/bottleneck.pyx b/src/python/gudhi/bottleneck.pyx
deleted file mode 100644
index af011e88..00000000
--- a/src/python/gudhi/bottleneck.pyx
+++ /dev/null
@@ -1,48 +0,0 @@
-# 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):       Vincent Rouvreau
-#
-# Copyright (C) 2016 Inria
-#
-# Modification(s):
-#   - YYYY/MM Author: Description of the modification
-
-from cython cimport numeric
-from libcpp.vector cimport vector
-from libcpp.utility cimport pair
-import os
-
-__author__ = "Vincent Rouvreau"
-__copyright__ = "Copyright (C) 2016 Inria"
-__license__ = "GPL v3"
-
-cdef extern from "Bottleneck_distance_interface.h" namespace "Gudhi::persistence_diagram":
-    double bottleneck(vector[pair[double, double]], vector[pair[double, double]], double)
-    double bottleneck(vector[pair[double, double]], vector[pair[double, double]])
-
-def bottleneck_distance(diagram_1, diagram_2, e=None):
-    """This function returns the point corresponding to a given vertex.
-
-    :param diagram_1: The first diagram.
-    :type diagram_1: vector[pair[double, double]]
-    :param diagram_2: The second diagram.
-    :type diagram_2: vector[pair[double, double]]
-    :param e: If `e` is 0, this uses an expensive algorithm to compute the
-        exact distance.
-        If `e` is not 0, it asks for an additive `e`-approximation, and
-        currently also allows a small multiplicative error (the last 2 or 3
-        bits of the mantissa may be wrong). This version of the algorithm takes
-        advantage of the limited precision of `double` and is usually a lot
-        faster to compute, whatever the value of `e`.
-
-        Thus, by default, `e` is the smallest positive double.
-    :type e: float
-    :rtype: float
-    :returns: the bottleneck distance.
-    """
-    if e is None:
-        # Default value is the smallest double value (not 0, 0 is for exact version)
-        return bottleneck(diagram_1, diagram_2)
-    else:
-        # Can be 0 for exact version
-        return bottleneck(diagram_1, diagram_2, e)
diff --git a/src/python/gudhi/hera.cc b/src/python/gudhi/hera.cc
index 5aec1806..ea80a9a8 100644
--- a/src/python/gudhi/hera.cc
+++ b/src/python/gudhi/hera.cc
@@ -8,39 +8,9 @@
  *      - YYYY/MM Author: Description of the modification
  */
 
-#include <pybind11/pybind11.h>
-#include <pybind11/numpy.h>
-
-#include <boost/range/counting_range.hpp>
-#include <boost/range/adaptor/transformed.hpp>
-
 #include <wasserstein.h> // Hera
 
-#include <utility>
-
-namespace py = pybind11;
-typedef py::array_t<double> Dgm;
-
-// Get m[i,0] and m[i,1] as a pair
-static auto pairify(void* p, ssize_t h, ssize_t w) {
-  return [=](ssize_t i){
-    char* birth = (char*)p + i * h;
-    char* death = birth + w;
-    return std::make_pair(*(double*)birth, *(double*)death);
-  };
-}
-
-inline auto numpy_to_range_of_pairs(py::array_t<double> dgm) {
-  py::buffer_info buf = dgm.request();
-  // shape (n,2) or (0) for empty
-  if((buf.ndim!=2 || buf.shape[1]!=2) && (buf.ndim!=1 || buf.shape[0]!=0))
-    throw std::runtime_error("Diagram must be an array of size n x 2");
-  // In the case of shape (0), avoid reading non-existing strides[1] even if we won't use it.
-  ssize_t stride1 = buf.ndim == 2 ? buf.strides[1] : 0;
-  auto cnt = boost::counting_range<ssize_t>(0, buf.shape[0]);
-  return boost::adaptors::transform(cnt, pairify(buf.ptr, buf.strides[0], stride1));
-  // Be careful that the returned range cannot contain references to dead temporaries.
-}
+#include <pybind11_diagram_utils.h>
 
 double wasserstein_distance(
     Dgm d1, Dgm d2,
diff --git a/src/python/include/pybind11_diagram_utils.h b/src/python/include/pybind11_diagram_utils.h
new file mode 100644
index 00000000..d9627258
--- /dev/null
+++ b/src/python/include/pybind11_diagram_utils.h
@@ -0,0 +1,39 @@
+/*    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):       Marc Glisse
+ *
+ *    Copyright (C) 2020 Inria
+ *
+ *    Modification(s):
+ *      - YYYY/MM Author: Description of the modification
+ */
+
+#include <pybind11/pybind11.h>
+#include <pybind11/numpy.h>
+
+#include <boost/range/counting_range.hpp>
+#include <boost/range/adaptor/transformed.hpp>
+
+namespace py = pybind11;
+typedef py::array_t<double> Dgm;
+
+// Get m[i,0] and m[i,1] as a pair
+static auto pairify(void* p, ssize_t h, ssize_t w) {
+  return [=](ssize_t i){
+    char* birth = (char*)p + i * h;
+    char* death = birth + w;
+    return std::make_pair(*(double*)birth, *(double*)death);
+  };
+}
+
+inline auto numpy_to_range_of_pairs(py::array_t<double> dgm) {
+  py::buffer_info buf = dgm.request();
+  // shape (n,2) or (0) for empty
+  if((buf.ndim!=2 || buf.shape[1]!=2) && (buf.ndim!=1 || buf.shape[0]!=0))
+    throw std::runtime_error("Diagram must be an array of size n x 2");
+  // In the case of shape (0), avoid reading non-existing strides[1] even if we won't use it.
+  ssize_t stride1 = buf.ndim == 2 ? buf.strides[1] : 0;
+  auto cnt = boost::counting_range<ssize_t>(0, buf.shape[0]);
+  return boost::adaptors::transform(cnt, pairify(buf.ptr, buf.strides[0], stride1));
+  // Be careful that the returned range cannot contain references to dead temporaries.
+}
diff --git a/src/python/setup.py.in b/src/python/setup.py.in
index f968bd59..852da910 100644
--- a/src/python/setup.py.in
+++ b/src/python/setup.py.in
@@ -18,7 +18,8 @@ __author__ = "Vincent Rouvreau"
 __copyright__ = "Copyright (C) 2016  Inria"
 __license__ = "MIT"
 
-modules = [@GUDHI_PYTHON_MODULES_TO_COMPILE@]
+cython_modules = [@GUDHI_CYTHON_MODULES@]
+pybind11_modules = [@GUDHI_PYBIND11_MODULES@]
 
 source_dir='@CMAKE_CURRENT_SOURCE_DIR@/gudhi/'
 extra_compile_args=[@GUDHI_PYTHON_EXTRA_COMPILE_ARGS@]
@@ -30,7 +31,7 @@ runtime_library_dirs=[@GUDHI_PYTHON_RUNTIME_LIBRARY_DIRS@]
 
 # Create ext_modules list from module list
 ext_modules = []
-for module in modules:
+for module in cython_modules:
     ext_modules.append(Extension(
         'gudhi.' + module,
         sources = [source_dir + module + '.pyx',],
@@ -55,6 +56,20 @@ ext_modules.append(Extension(
     extra_compile_args=extra_compile_args + [@GUDHI_PYBIND11_EXTRA_COMPILE_ARGS@],
     ))
 
+if "bottleneck" in pybind11_modules:
+    ext_modules.append(Extension(
+        'gudhi.bottleneck',
+        sources = [source_dir + 'bottleneck.cc'],
+        language = 'c++',
+        include_dirs = include_dirs +
+            [pybind11.get_include(False), pybind11.get_include(True)],
+        extra_compile_args=extra_compile_args + [@GUDHI_PYBIND11_EXTRA_COMPILE_ARGS@],
+        extra_link_args=extra_link_args,
+        libraries=libraries,
+        library_dirs=library_dirs,
+        runtime_library_dirs=runtime_library_dirs,
+        ))
+
 setup(
     name = 'gudhi',
     packages=find_packages(), # find_namespace_packages(include=["gudhi*"])
-- 
cgit v1.2.3


From d61bfd349274456f8d7e0ccd64839a2d84eea0a0 Mon Sep 17 00:00:00 2001
From: Marc Glisse <marc.glisse@inria.fr>
Date: Thu, 7 May 2020 08:40:55 +0200
Subject: doc

---
 src/python/gudhi/bottleneck.cc | 5 ++---
 1 file changed, 2 insertions(+), 3 deletions(-)

diff --git a/src/python/gudhi/bottleneck.cc b/src/python/gudhi/bottleneck.cc
index 577e5e0b..732cb9a8 100644
--- a/src/python/gudhi/bottleneck.cc
+++ b/src/python/gudhi/bottleneck.cc
@@ -32,9 +32,9 @@ PYBIND11_MODULE(bottleneck, m) {
     This function returns the point corresponding to a given vertex.
 
     :param diagram_1: The first diagram.
-    :type diagram_1: vector[pair[double, double]]
+    :type diagram_1: numpy array of shape (m,2)
     :param diagram_2: The second diagram.
-    :type diagram_2: vector[pair[double, double]]
+    :type diagram_2: numpy array of shape (n,2)
     :param e: If `e` is 0, this uses an expensive algorithm to compute the
         exact distance.
         If `e` is not 0, it asks for an additive `e`-approximation, and
@@ -42,7 +42,6 @@ PYBIND11_MODULE(bottleneck, m) {
         bits of the mantissa may be wrong). This version of the algorithm takes
         advantage of the limited precision of `double` and is usually a lot
         faster to compute, whatever the value of `e`.
-
         Thus, by default, `e` is the smallest positive double.
     :type e: float
     :rtype: float
-- 
cgit v1.2.3


From acc76eb90b8cfe3f8cbb8d30f101c7f879ab61c4 Mon Sep 17 00:00:00 2001
From: Marc Glisse <marc.glisse@inria.fr>
Date: Thu, 7 May 2020 20:10:46 +0200
Subject: Warn for initialize_filtration

---
 src/python/gudhi/simplex_tree.pyx | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/src/python/gudhi/simplex_tree.pyx b/src/python/gudhi/simplex_tree.pyx
index 55115cca..b23885b4 100644
--- a/src/python/gudhi/simplex_tree.pyx
+++ b/src/python/gudhi/simplex_tree.pyx
@@ -101,6 +101,8 @@ cdef class SimplexTree:
 
         .. deprecated:: 3.2.0
         """
+        import warnings
+        warnings.warn("Since Gudhi 3.2, calling SimplexTree.initialize_filtration is unnecessary.", DeprecationWarning)
         self.get_ptr().initialize_filtration()
 
     def num_vertices(self):
-- 
cgit v1.2.3


From 778c0af7dea0c103db85986fe2e2eb5fddd7588f Mon Sep 17 00:00:00 2001
From: Marc Glisse <marc.glisse@inria.fr>
Date: Fri, 8 May 2020 10:14:50 +0200
Subject: Loop on pybind11 modules

---
 src/python/setup.py.in | 22 +++++++---------------
 1 file changed, 7 insertions(+), 15 deletions(-)

diff --git a/src/python/setup.py.in b/src/python/setup.py.in
index 852da910..b9f4e3f0 100644
--- a/src/python/setup.py.in
+++ b/src/python/setup.py.in
@@ -46,23 +46,15 @@ for module in cython_modules:
 
 ext_modules = cythonize(ext_modules)
 
-ext_modules.append(Extension(
-    'gudhi.hera',
-    sources = [source_dir + 'hera.cc'],
-    language = 'c++',
-    include_dirs = include_dirs +
-        ['@HERA_WASSERSTEIN_INCLUDE_DIR@',
-        pybind11.get_include(False), pybind11.get_include(True)],
-    extra_compile_args=extra_compile_args + [@GUDHI_PYBIND11_EXTRA_COMPILE_ARGS@],
-    ))
-
-if "bottleneck" in pybind11_modules:
+for module in pybind11_modules:
+    my_include_dirs = include_dirs + [pybind11.get_include(False), pybind11.get_include(True)]
+    if module == 'hera':
+        my_include_dirs = ['@HERA_WASSERSTEIN_INCLUDE_DIR@'] + my_include_dirs
     ext_modules.append(Extension(
-        'gudhi.bottleneck',
-        sources = [source_dir + 'bottleneck.cc'],
+        'gudhi.' + module,
+        sources = [source_dir + module + '.cc'],
         language = 'c++',
-        include_dirs = include_dirs +
-            [pybind11.get_include(False), pybind11.get_include(True)],
+        include_dirs = my_include_dirs,
         extra_compile_args=extra_compile_args + [@GUDHI_PYBIND11_EXTRA_COMPILE_ARGS@],
         extra_link_args=extra_link_args,
         libraries=libraries,
-- 
cgit v1.2.3