added call methods + other fixes

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,:]