merging master (including wasserstein distance) into wbary ; modified CMakeList to solve conflicts

20 files changed, 413 insertions, 93 deletions

diff --git a/src/Simplex_tree/example/simple_simplex_tree.cpp b/src/Simplex_tree/example/simple_simplex_tree.cpp
index 4353939f..47ea7e36 100644
--- a/src/Simplex_tree/example/simple_simplex_tree.cpp
+++ b/src/Simplex_tree/example/simple_simplex_tree.cpp
@@ -166,10 +166,19 @@ int main(int argc, char* const argv[]) {
   // ++ GENERAL VARIABLE SET
 
   std::cout << "********************************************************************\n";
-  // Display the Simplex_tree - Can not be done in the middle of 2 inserts
   std::cout << "* The complex contains " << simplexTree.num_simplices() << " simplices\n";
   std::cout << "   - dimension " << simplexTree.dimension() << "\n";
-  std::cout << "* Iterator on Simplices in the filtration, with [filtration value]:\n";
+  std::cout << "* Iterator on simplices, with [filtration value]:\n";
+  for (Simplex_tree::Simplex_handle f_simplex : simplexTree.complex_simplex_range()) {
+    std::cout << "   "
+              << "[" << simplexTree.filtration(f_simplex) << "] ";
+    for (auto vertex : simplexTree.simplex_vertex_range(f_simplex)) std::cout << "(" << vertex << ")";
+    std::cout << std::endl;
+  }
+
+  std::cout << "********************************************************************\n";
+  // Can not be done in the middle of 2 inserts
+  std::cout << "* Iterator on simplices sorted by filtration values, with [filtration value]:\n";
   for (auto f_simplex : simplexTree.filtration_simplex_range()) {
     std::cout << "   "
               << "[" << simplexTree.filtration(f_simplex) << "] ";
diff --git a/src/python/CMakeLists.txt b/src/python/CMakeLists.txt
index 90c399b4..9fa7b129 100644
--- a/src/python/CMakeLists.txt
+++ b/src/python/CMakeLists.txt
@@ -57,6 +57,7 @@ if(PYTHONINTERP_FOUND)
     set(GUDHI_PYTHON_MODULES_EXTRA "${GUDHI_PYTHON_MODULES_EXTRA}'representations', ")
     set(GUDHI_PYTHON_MODULES_EXTRA "${GUDHI_PYTHON_MODULES_EXTRA}'wasserstein', ")
     set(GUDHI_PYTHON_MODULES_EXTRA "${GUDHI_PYTHON_MODULES_EXTRA}'barycenter', ")
+    set(GUDHI_PYTHON_MODULES_EXTRA "${GUDHI_PYTHON_MODULES_EXTRA}'point_cloud', ")
 
     add_gudhi_debug_info("Python version ${PYTHON_VERSION_STRING}")
     add_gudhi_debug_info("Cython version ${CYTHON_VERSION}")
@@ -228,6 +229,7 @@ if(PYTHONINTERP_FOUND)
     file(COPY "gudhi/representations" DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/gudhi/")
     file(COPY "gudhi/wasserstein.py" DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/gudhi")
     file(COPY "gudhi/barycenter.py" DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/gudhi")
+    file(COPY "gudhi/point_cloud" DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/gudhi")
 
     add_custom_command(
         OUTPUT gudhi.so
@@ -407,6 +409,9 @@ if(PYTHONINTERP_FOUND)
       add_gudhi_py_test(test_representations)
     endif()
 
+    # Time Delay
+    add_gudhi_py_test(test_time_delay)
+
     # Documentation generation is available through sphinx - requires all modules
     if(SPHINX_PATH)
       if(MATPLOTLIB_FOUND)
diff --git a/src/python/doc/point_cloud.rst b/src/python/doc/point_cloud.rst
index d668428a..c0d4b303 100644
--- a/src/python/doc/point_cloud.rst
+++ b/src/python/doc/point_cloud.rst
@@ -20,3 +20,11 @@ Subsampling
    :members:
    :special-members:
    :show-inheritance:
+
+TimeDelayEmbedding
+------------------
+
+.. autoclass:: gudhi.point_cloud.timedelay.TimeDelayEmbedding
+   :members:
+   :special-members: __call__
+
diff --git a/src/python/doc/wasserstein_distance_user.rst b/src/python/doc/wasserstein_distance_user.rst
index 94b454e2..a9b21fa5 100644
--- a/src/python/doc/wasserstein_distance_user.rst
+++ b/src/python/doc/wasserstein_distance_user.rst
@@ -36,10 +36,10 @@ Note that persistence diagrams must be submitted as (n x 2) numpy arrays and mus
     import gudhi.wasserstein
     import numpy as np
 
-    diag1 = np.array([[2.7, 3.7],[9.6, 14.],[34.2, 34.974]])
-    diag2 = np.array([[2.8, 4.45],[9.5, 14.1]])
+    dgm1 = np.array([[2.7, 3.7],[9.6, 14.],[34.2, 34.974]])
+    dgm2 = np.array([[2.8, 4.45],[9.5, 14.1]])
 
-    message = "Wasserstein distance value = " + '%.2f' % gudhi.wasserstein.wasserstein_distance(diag1, diag2, order=1., internal_p=2.)
+    message = "Wasserstein distance value = " + '%.2f' % gudhi.wasserstein.wasserstein_distance(dgm1, dgm2, order=1., internal_p=2.)
     print(message)
 
 The output is:
@@ -47,3 +47,40 @@ The output is:
 .. testoutput::
 
     Wasserstein distance value = 1.45
+
+We can also have access to the optimal matching by letting `matching=True`. 
+It is encoded as a list of indices (i,j), meaning that the i-th point in X
+is mapped to the j-th point in Y. 
+An index of -1 represents the diagonal.
+
+.. testcode::
+
+    import gudhi.wasserstein
+    import numpy as np
+
+    dgm1 = np.array([[2.7, 3.7],[9.6, 14.],[34.2, 34.974]])
+    dgm2 = np.array([[2.8, 4.45], [5, 6], [9.5, 14.1]])
+    cost, matchings = gudhi.wasserstein.wasserstein_distance(dgm1, dgm2, matching=True, order=1, internal_p=2)
+
+    message_cost = "Wasserstein distance value = %.2f" %cost
+    print(message_cost)
+    dgm1_to_diagonal = matchings[matchings[:,1] == -1, 0]
+    dgm2_to_diagonal = matchings[matchings[:,0] == -1, 1]
+    off_diagonal_match = np.delete(matchings, np.where(matchings == -1)[0], axis=0)
+
+    for i,j in off_diagonal_match:
+        print("point %s in dgm1 is matched to point %s in dgm2" %(i,j))
+    for i in dgm1_to_diagonal:
+        print("point %s in dgm1 is matched to the diagonal" %i)
+    for j in dgm2_to_diagonal:
+        print("point %s in dgm2 is matched to the diagonal" %j)
+
+The output is:
+
+.. testoutput::
+    
+    Wasserstein distance value = 2.15
+    point 0 in dgm1 is matched to point 0 in dgm2
+    point 1 in dgm1 is matched to point 2 in dgm2
+    point 2 in dgm1 is matched to the diagonal
+    point 1 in dgm2 is matched to the diagonal
diff --git a/src/python/example/alpha_complex_from_points_example.py b/src/python/example/alpha_complex_from_points_example.py
index 844d7a82..73faf17c 100755
--- a/src/python/example/alpha_complex_from_points_example.py
+++ b/src/python/example/alpha_complex_from_points_example.py
@@ -46,8 +46,14 @@ if simplex_tree.find([4]):
 else:
     print("[4] Not found...")
 
+# Some insertions, simplex_tree needs to initialize filtrations
+simplex_tree.initialize_filtration()
+
 print("dimension=", simplex_tree.dimension())
-print("filtrations=", simplex_tree.get_filtration())
+print("filtrations=")
+for simplex_with_filtration in simplex_tree.get_filtration():
+    print("(%s, %.2f)" % tuple(simplex_with_filtration))
+
 print("star([0])=", simplex_tree.get_star([0]))
 print("coface([0], 1)=", simplex_tree.get_cofaces([0], 1))
 
diff --git a/src/python/example/rips_complex_from_points_example.py b/src/python/example/rips_complex_from_points_example.py
index 59d8a261..c05703c6 100755
--- a/src/python/example/rips_complex_from_points_example.py
+++ b/src/python/example/rips_complex_from_points_example.py
@@ -22,6 +22,9 @@ rips = gudhi.RipsComplex(points=[[0, 0], [1, 0], [0, 1], [1, 1]], max_edge_lengt
 
 simplex_tree = rips.create_simplex_tree(max_dimension=1)
 
-print("filtrations=", simplex_tree.get_filtration())
+print("filtrations=")
+for simplex_with_filtration in simplex_tree.get_filtration():
+    print("(%s, %.2f)" % tuple(simplex_with_filtration))
+
 print("star([0])=", simplex_tree.get_star([0]))
 print("coface([0], 1)=", simplex_tree.get_cofaces([0], 1))
diff --git a/src/python/example/simplex_tree_example.py b/src/python/example/simplex_tree_example.py
index 30de00da..34833899 100755
--- a/src/python/example/simplex_tree_example.py
+++ b/src/python/example/simplex_tree_example.py
@@ -38,8 +38,15 @@ else:
 
 print("dimension=", st.dimension())
 
+print("simplices=")
+for simplex_with_filtration in st.get_simplices():
+    print("(%s, %.2f)" % tuple(simplex_with_filtration))
+
 st.initialize_filtration()
-print("filtration=", st.get_filtration())
+print("filtration=")
+for simplex_with_filtration in st.get_filtration():
+    print("(%s, %.2f)" % tuple(simplex_with_filtration))
+
 print("filtration[1, 2]=", st.filtration([1, 2]))
 print("filtration[4, 2]=", st.filtration([4, 2]))
 
diff --git a/src/python/gudhi/point_cloud/__init__.py b/src/python/gudhi/point_cloud/__init__.py
new file mode 100644
index 00000000..e69de29b
--- /dev/null
+++ b/src/python/gudhi/point_cloud/__init__.py
diff --git a/src/python/gudhi/point_cloud/timedelay.py b/src/python/gudhi/point_cloud/timedelay.py
new file mode 100644
index 00000000..f01df442
--- /dev/null
+++ b/src/python/gudhi/point_cloud/timedelay.py
@@ -0,0 +1,95 @@
+# 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):       Martin Royer, Yuichi Ike, Masatoshi Takenouchi
+#
+# Copyright (C) 2020 Inria, Copyright (C) 2020 Fujitsu Laboratories Ltd.
+# Modification(s):
+#   - YYYY/MM Author: Description of the modification
+
+import numpy as np
+
+
+class TimeDelayEmbedding:
+    """Point cloud transformation class.
+    Embeds time-series data in the R^d according to [Takens' Embedding Theorem]
+    (https://en.wikipedia.org/wiki/Takens%27s_theorem) and obtains the 
+    coordinates of each point.
+
+    Parameters
+    ----------
+    dim : int, optional (default=3)
+        `d` of R^d to be embedded.
+    delay : int, optional (default=1)
+        Time-Delay embedding.
+    skip : int, optional (default=1)
+        How often to skip embedded points.
+
+    Example
+    -------
+
+    Given delay=3 and skip=2, a point cloud which is obtained by embedding
+    a scalar time-series into R^3 is as follows::
+
+        time-series = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
+        point cloud = [[1, 4, 7],
+                       [3, 6, 9]]
+
+    Given delay=1 and skip=1, a point cloud which is obtained by embedding
+    a 2D vector time-series data into R^4 is as follows::
+
+        time-series = [[0, 1], [2, 3], [4, 5], [6, 7], [8, 9]]
+        point cloud = [[0, 1, 2, 3],
+                       [2, 3, 4, 5],
+                       [4, 5, 6, 7],
+                       [6, 7, 8, 9]]
+    """
+
+    def __init__(self, dim=3, delay=1, skip=1):
+        self._dim = dim
+        self._delay = delay
+        self._skip = skip
+
+    def __call__(self, ts):
+        """Transform method for single time-series data.
+
+        Parameters
+        ----------
+        ts : Iterable[float] or Iterable[Iterable[float]]
+            A single time-series data, with scalar or vector values.
+
+        Returns
+        -------
+        point cloud : n x dim numpy arrays
+            Makes point cloud from a single time-series data.
+        """
+        return self._transform(np.array(ts))
+
+    def fit(self, ts, y=None):
+        return self
+    
+    def _transform(self, ts):
+        """Guts of transform method."""
+        if ts.ndim == 1:
+            repeat = self._dim
+        else:
+            assert self._dim % ts.shape[1] == 0
+            repeat = self._dim // ts.shape[1]
+        end = len(ts) - self._delay * (repeat - 1)
+        short = np.arange(0, end, self._skip)
+        vertical = np.arange(0, repeat * self._delay, self._delay)
+        return ts[np.add.outer(short, vertical)].reshape(len(short), -1)
+
+    def transform(self, ts):
+        """Transform method for multiple time-series data.
+
+        Parameters
+        ----------
+        ts : Iterable[Iterable[float]] or Iterable[Iterable[Iterable[float]]]
+            Multiple time-series data, with scalar or vector values.
+
+        Returns
+        -------
+        point clouds : list of n x dim numpy arrays
+            Makes point cloud from each time-series data.
+        """
+        return [self._transform(np.array(s)) for s in ts]
diff --git a/src/python/gudhi/simplex_tree.pxd b/src/python/gudhi/simplex_tree.pxd
index 96d14079..82f155de 100644
--- a/src/python/gudhi/simplex_tree.pxd
+++ b/src/python/gudhi/simplex_tree.pxd
@@ -21,6 +21,22 @@ cdef extern from "Simplex_tree_interface.h" namespace "Gudhi":
     cdef cppclass Simplex_tree_options_full_featured:
         pass
 
+    cdef cppclass Simplex_tree_simplex_handle "Gudhi::Simplex_tree_interface<Gudhi::Simplex_tree_options_full_featured>::Simplex_handle":
+        pass
+
+    cdef cppclass Simplex_tree_simplices_iterator "Gudhi::Simplex_tree_interface<Gudhi::Simplex_tree_options_full_featured>::Complex_simplex_iterator":
+        Simplex_tree_simplices_iterator()
+        Simplex_tree_simplex_handle& operator*()
+        Simplex_tree_simplices_iterator operator++()
+        bint operator!=(Simplex_tree_simplices_iterator)
+
+    cdef cppclass Simplex_tree_skeleton_iterator "Gudhi::Simplex_tree_interface<Gudhi::Simplex_tree_options_full_featured>::Skeleton_simplex_iterator":
+        Simplex_tree_skeleton_iterator()
+        Simplex_tree_simplex_handle& operator*()
+        Simplex_tree_skeleton_iterator operator++()
+        bint operator!=(Simplex_tree_skeleton_iterator)
+
+
     cdef cppclass Simplex_tree_interface_full_featured "Gudhi::Simplex_tree_interface<Gudhi::Simplex_tree_options_full_featured>":
         Simplex_tree()
         double simplex_filtration(vector[int] simplex)
@@ -34,8 +50,6 @@ cdef extern from "Simplex_tree_interface.h" namespace "Gudhi":
         bool find_simplex(vector[int] simplex)
         bool insert_simplex_and_subfaces(vector[int] simplex,
                                          double filtration)
-        vector[pair[vector[int], double]] get_filtration()
-        vector[pair[vector[int], double]] get_skeleton(int dimension)
         vector[pair[vector[int], double]] get_star(vector[int] simplex)
         vector[pair[vector[int], double]] get_cofaces(vector[int] simplex,
                                                           int dimension)
@@ -43,6 +57,14 @@ cdef extern from "Simplex_tree_interface.h" namespace "Gudhi":
         void remove_maximal_simplex(vector[int] simplex)
         bool prune_above_filtration(double filtration)
         bool make_filtration_non_decreasing()
+        # Iterators over Simplex tree
+        pair[vector[int], double] get_simplex_and_filtration(Simplex_tree_simplex_handle f_simplex)
+        Simplex_tree_simplices_iterator get_simplices_iterator_begin()
+        Simplex_tree_simplices_iterator get_simplices_iterator_end()
+        vector[Simplex_tree_simplex_handle].const_iterator get_filtration_iterator_begin()
+        vector[Simplex_tree_simplex_handle].const_iterator get_filtration_iterator_end()
+        Simplex_tree_skeleton_iterator get_skeleton_iterator_begin(int dimension)
+        Simplex_tree_skeleton_iterator get_skeleton_iterator_end(int dimension)
 
 cdef extern from "Persistent_cohomology_interface.h" namespace "Gudhi":
     cdef cppclass Simplex_tree_persistence_interface "Gudhi::Persistent_cohomology_interface<Gudhi::Simplex_tree<Gudhi::Simplex_tree_options_full_featured>>":
diff --git a/src/python/gudhi/simplex_tree.pyx b/src/python/gudhi/simplex_tree.pyx
index b18627c4..c01cc905 100644
--- a/src/python/gudhi/simplex_tree.pyx
+++ b/src/python/gudhi/simplex_tree.pyx
@@ -7,6 +7,7 @@
 # Modification(s):
 #   - YYYY/MM Author: Description of the modification
 
+from cython.operator import dereference, preincrement
 from libc.stdint cimport intptr_t
 from numpy import array as np_array
 cimport simplex_tree
@@ -207,22 +208,34 @@ cdef class SimplexTree:
         return self.get_ptr().insert_simplex_and_subfaces(csimplex,
                                                         <double>filtration)
 
-    def get_filtration(self):
-        """This function returns a list of all simplices with their given
+    def get_simplices(self):
+        """This function returns a generator with simplices and their given
         filtration values.
 
+        :returns:  The simplices.
+        :rtype:  generator with tuples(simplex, filtration)
+        """
+        cdef Simplex_tree_simplices_iterator it = self.get_ptr().get_simplices_iterator_begin()
+        cdef Simplex_tree_simplices_iterator end = self.get_ptr().get_simplices_iterator_end()
+        cdef Simplex_tree_simplex_handle sh = dereference(it)
+
+        while it != end:
+            yield self.get_ptr().get_simplex_and_filtration(dereference(it))
+            preincrement(it)
+
+    def get_filtration(self):
+        """This function returns a generator with simplices and their given
+        filtration values sorted by increasing filtration values.
+
         :returns:  The simplices sorted by increasing filtration values.
-        :rtype:  list of tuples(simplex, filtration)
+        :rtype:  generator with tuples(simplex, filtration)
         """
-        cdef vector[pair[vector[int], double]] filtration \
-            = self.get_ptr().get_filtration()
-        ct = []
-        for filtered_complex in filtration:
-            v = []
-            for vertex in filtered_complex.first:
-                v.append(vertex)
-            ct.append((v, filtered_complex.second))
-        return ct
+        cdef vector[Simplex_tree_simplex_handle].const_iterator it = self.get_ptr().get_filtration_iterator_begin()
+        cdef vector[Simplex_tree_simplex_handle].const_iterator end = self.get_ptr().get_filtration_iterator_end()
+
+        while it != end:
+            yield self.get_ptr().get_simplex_and_filtration(dereference(it))
+            preincrement(it)
 
     def get_skeleton(self, dimension):
         """This function returns the (simplices of the) skeleton of a maximum
@@ -233,15 +246,12 @@ cdef class SimplexTree:
         :returns:  The (simplices of the) skeleton of a maximum dimension.
         :rtype:  list of tuples(simplex, filtration)
         """
-        cdef vector[pair[vector[int], double]] skeleton \
-            = self.get_ptr().get_skeleton(<int>dimension)
-        ct = []
-        for filtered_simplex in skeleton:
-            v = []
-            for vertex in filtered_simplex.first:
-                v.append(vertex)
-            ct.append((v, filtered_simplex.second))
-        return ct
+        cdef Simplex_tree_skeleton_iterator it = self.get_ptr().get_skeleton_iterator_begin(dimension)
+        cdef Simplex_tree_skeleton_iterator end = self.get_ptr().get_skeleton_iterator_end(dimension)
+
+        while it != end:
+            yield self.get_ptr().get_simplex_and_filtration(dereference(it))
+            preincrement(it)
 
     def get_star(self, simplex):
         """This function returns the star of a given N-simplex.
diff --git a/src/python/gudhi/wasserstein.py b/src/python/gudhi/wasserstein.py
index 13102094..3dd993f9 100644
--- a/src/python/gudhi/wasserstein.py
+++ b/src/python/gudhi/wasserstein.py
@@ -30,8 +30,10 @@ def _build_dist_matrix(X, Y, order=2., internal_p=2.):
     :param order: exponent for the Wasserstein metric.
     :param internal_p: Ground metric (i.e. norm L^p).
     :returns: (n+1) x (m+1) np.array encoding the cost matrix C. 
-                For 1 <= i <= n, 1 <= j <= m, C[i,j] encodes the distance between X[i] and Y[j], while C[i, m+1] (resp. C[n+1, j]) encodes the distance (to the p) between X[i] (resp Y[j]) and its orthogonal proj onto the diagonal.
-                note also that C[n+1, m+1] = 0  (it costs nothing to move from the diagonal to the diagonal).
+                For 0 <= i < n, 0 <= j < m, C[i,j] encodes the distance between X[i] and Y[j], 
+                while C[i, m] (resp. C[n, j]) encodes the distance (to the p) between X[i] (resp Y[j]) 
+                and its orthogonal projection onto the diagonal.
+                note also that C[n, m] = 0  (it costs nothing to move from the diagonal to the diagonal).
     '''
     Xdiag = _proj_on_diag(X)
     Ydiag = _proj_on_diag(Y)
@@ -62,14 +64,20 @@ def _perstot(X, order, internal_p):
     return (np.sum(np.linalg.norm(X - Xdiag, ord=internal_p, axis=1)**order))**(1./order)
 
 
-def wasserstein_distance(X, Y, order=2., internal_p=2.):
+def wasserstein_distance(X, Y, matching=False, order=2., internal_p=2.):
     '''
-    :param X: (n x 2) numpy.array encoding the (finite points of the) first diagram. Must not contain essential points (i.e. with infinite coordinate).
+    :param X: (n x 2) numpy.array encoding the (finite points of the) first diagram. Must not contain essential points 
+                (i.e. with infinite coordinate).
     :param Y: (m x 2) numpy.array encoding the second diagram.
+    :param matching: if True, computes and returns the optimal matching between X and Y, encoded as
+                     a (n x 2) np.array  [...[i,j]...], meaning the i-th point in X is matched to
+                     the j-th point in Y, with the convention (-1) represents the diagonal.
     :param order: exponent for Wasserstein; Default value is 2.
-    :param internal_p: Ground metric on the (upper-half) plane (i.e. norm L^p in R^2); Default value is 2 (Euclidean norm).
-    :returns: the Wasserstein distance of order q (1 <= q < infinity) between persistence diagrams with respect to the internal_p-norm as ground metric.
-    :rtype: float
+    :param internal_p: Ground metric on the (upper-half) plane (i.e. norm L^p in R^2); 
+                       Default value is 2 (Euclidean norm).
+    :returns: the Wasserstein distance of order q (1 <= q < infinity) between persistence diagrams with 
+              respect to the internal_p-norm as ground metric.
+              If matching is set to True, also returns the optimal matching between X and Y.
     '''
     n = len(X)
     m = len(Y)
@@ -77,21 +85,40 @@ def wasserstein_distance(X, Y, order=2., internal_p=2.):
     # handle empty diagrams
     if X.size == 0:
         if Y.size == 0:
-            return 0.
+            if not matching:
+                return 0.
+            else:
+                return 0., np.array([])
         else:
-            return _perstot(Y, order, internal_p)
+            if not matching:
+                return _perstot(Y, order, internal_p)
+            else:
+                return _perstot(Y, order, internal_p), np.array([[-1, j] for j in range(m)])
     elif Y.size == 0:
-        return _perstot(X, order, internal_p)
+        if not matching:
+            return _perstot(X, order, internal_p)
+        else:
+            return _perstot(X, order, internal_p), np.array([[i, -1] for i in range(n)])
 
     M = _build_dist_matrix(X, Y, order=order, internal_p=internal_p)
-    a = np.full(n+1, 1. / (n + m) )  # weight vector of the input diagram. Uniform here.
-    a[-1] = a[-1] * m                # normalized so that we have a probability measure, required by POT
-    b = np.full(m+1, 1. / (n + m) )  # weight vector of the input diagram. Uniform here.
-    b[-1] = b[-1] * n                # so that we have a probability measure, required by POT
+    a = np.ones(n+1) # weight vector of the input diagram. Uniform here.
+    a[-1] = m
+    b = np.ones(m+1) # weight vector of the input diagram. Uniform here.
+    b[-1] = n
+
+    if matching:
+        P = ot.emd(a=a,b=b,M=M, numItermax=2000000)
+        ot_cost = np.sum(np.multiply(P,M))
+        P[-1, -1] = 0  # Remove matching corresponding to the diagonal
+        match = np.argwhere(P)
+        # Now we turn to -1 points encoding the diagonal
+        match[:,0][match[:,0] >= n] = -1
+        match[:,1][match[:,1] >= m] = -1
+        return ot_cost ** (1./order) , match
 
     # Comptuation of the otcost using the ot.emd2 library.
     # Note: it is the Wasserstein distance to the power q.
     # The default numItermax=100000 is not sufficient for some examples with 5000 points, what is a good value?
-    ot_cost = (n+m) * ot.emd2(a, b, M, numItermax=2000000)
+    ot_cost = ot.emd2(a, b, M, numItermax=2000000)
 
     return ot_cost ** (1./order)
diff --git a/src/python/include/Simplex_tree_interface.h b/src/python/include/Simplex_tree_interface.h
index 06f31341..4a7062d6 100644
--- a/src/python/include/Simplex_tree_interface.h
+++ b/src/python/include/Simplex_tree_interface.h
@@ -33,7 +33,10 @@ class Simplex_tree_interface : public Simplex_tree<SimplexTreeOptions> {
   using Simplex_handle = typename Base::Simplex_handle;
   using Insertion_result = typename std::pair<Simplex_handle, bool>;
   using Simplex = std::vector<Vertex_handle>;
-  using Filtered_simplices = std::vector<std::pair<Simplex, Filtration_value>>;
+  using Simplex_and_filtration = std::pair<Simplex, Filtration_value>;
+  using Filtered_simplices = std::vector<Simplex_and_filtration>;
+  using Skeleton_simplex_iterator = typename Base::Skeleton_simplex_iterator;
+  using Complex_simplex_iterator = typename Base::Complex_simplex_iterator;
 
  public:
   bool find_simplex(const Simplex& vh) {
@@ -82,29 +85,12 @@ class Simplex_tree_interface : public Simplex_tree<SimplexTreeOptions> {
     Base::initialize_filtration();
   }
 
-  Filtered_simplices get_filtration() {
-    Base::initialize_filtration();
-    Filtered_simplices filtrations;
-    for (auto f_simplex : Base::filtration_simplex_range()) {
-      Simplex simplex;
-      for (auto vertex : Base::simplex_vertex_range(f_simplex)) {
-        simplex.insert(simplex.begin(), vertex);
-      }
-      filtrations.push_back(std::make_pair(simplex, Base::filtration(f_simplex)));
-    }
-    return filtrations;
-  }
-
-  Filtered_simplices get_skeleton(int dimension) {
-    Filtered_simplices skeletons;
-    for (auto f_simplex : Base::skeleton_simplex_range(dimension)) {
-      Simplex simplex;
-      for (auto vertex : Base::simplex_vertex_range(f_simplex)) {
-        simplex.insert(simplex.begin(), vertex);
-      }
-      skeletons.push_back(std::make_pair(simplex, Base::filtration(f_simplex)));
+  Simplex_and_filtration get_simplex_and_filtration(Simplex_handle f_simplex) {
+    Simplex simplex;
+    for (auto vertex : Base::simplex_vertex_range(f_simplex)) {
+      simplex.insert(simplex.begin(), vertex);
     }
-    return skeletons;
+    return std::make_pair(std::move(simplex), Base::filtration(f_simplex));
   }
 
   Filtered_simplices get_star(const Simplex& simplex) {
@@ -135,6 +121,38 @@ class Simplex_tree_interface : public Simplex_tree<SimplexTreeOptions> {
     Base::initialize_filtration();
     pcoh = new Gudhi::Persistent_cohomology_interface<Base>(*this);
   }
+
+  // Iterator over the simplex tree
+  Complex_simplex_iterator get_simplices_iterator_begin() {
+    // this specific case works because the range is just a pair of iterators - won't work if range was a vector
+    return Base::complex_simplex_range().begin();
+  }
+
+  Complex_simplex_iterator get_simplices_iterator_end() {
+    // this specific case works because the range is just a pair of iterators - won't work if range was a vector
+    return Base::complex_simplex_range().end();
+  }
+
+  typename std::vector<Simplex_handle>::const_iterator get_filtration_iterator_begin() {
+    // Base::initialize_filtration(); already performed in filtration_simplex_range
+    // this specific case works because the range is just a pair of iterators - won't work if range was a vector
+    return Base::filtration_simplex_range().begin();
+  }
+
+  typename std::vector<Simplex_handle>::const_iterator get_filtration_iterator_end() {
+    // this specific case works because the range is just a pair of iterators - won't work if range was a vector
+    return Base::filtration_simplex_range().end();
+  }
+
+  Skeleton_simplex_iterator get_skeleton_iterator_begin(int dimension) {
+    // this specific case works because the range is just a pair of iterators - won't work if range was a vector
+    return Base::skeleton_simplex_range(dimension).begin();
+  }
+
+  Skeleton_simplex_iterator get_skeleton_iterator_end(int dimension) {
+    // this specific case works because the range is just a pair of iterators - won't work if range was a vector
+    return Base::skeleton_simplex_range(dimension).end();
+  }
 };
 
 }  // namespace Gudhi
diff --git a/src/python/test/test_alpha_complex.py b/src/python/test/test_alpha_complex.py
index 3761fe16..77121302 100755
--- a/src/python/test/test_alpha_complex.py
+++ b/src/python/test/test_alpha_complex.py
@@ -40,7 +40,7 @@ def test_infinite_alpha():
     assert simplex_tree.num_simplices() == 11
     assert simplex_tree.num_vertices() == 4
 
-    assert simplex_tree.get_filtration() == [
+    assert list(simplex_tree.get_filtration()) == [
         ([0], 0.0),
         ([1], 0.0),
         ([2], 0.0),
@@ -53,6 +53,7 @@ def test_infinite_alpha():
         ([0, 1, 2], 0.5),
         ([1, 2, 3], 0.5),
     ]
+
     assert simplex_tree.get_star([0]) == [
         ([0], 0.0),
         ([0, 1], 0.25),
@@ -105,7 +106,7 @@ def test_filtered_alpha():
     else:
         assert False
 
-    assert simplex_tree.get_filtration() == [
+    assert list(simplex_tree.get_filtration()) == [
         ([0], 0.0),
         ([1], 0.0),
         ([2], 0.0),
diff --git a/src/python/test/test_euclidean_witness_complex.py b/src/python/test/test_euclidean_witness_complex.py
index c18d2484..f3664d39 100755
--- a/src/python/test/test_euclidean_witness_complex.py
+++ b/src/python/test/test_euclidean_witness_complex.py
@@ -40,7 +40,7 @@ def test_witness_complex():
     assert landmarks[1] == euclidean_witness_complex.get_point(1)
     assert landmarks[2] == euclidean_witness_complex.get_point(2)
 
-    assert simplex_tree.get_filtration() == [
+    assert list(simplex_tree.get_filtration()) == [
         ([0], 0.0),
         ([1], 0.0),
         ([0, 1], 0.0),
@@ -78,13 +78,13 @@ def test_strong_witness_complex():
     assert landmarks[1] == euclidean_strong_witness_complex.get_point(1)
     assert landmarks[2] == euclidean_strong_witness_complex.get_point(2)
 
-    assert simplex_tree.get_filtration() == [([0], 0.0), ([1], 0.0), ([2], 0.0)]
+    assert list(simplex_tree.get_filtration()) == [([0], 0.0), ([1], 0.0), ([2], 0.0)]
 
     simplex_tree = euclidean_strong_witness_complex.create_simplex_tree(
         max_alpha_square=100.0
     )
 
-    assert simplex_tree.get_filtration() == [
+    assert list(simplex_tree.get_filtration()) == [
         ([0], 0.0),
         ([1], 0.0),
         ([2], 0.0),
diff --git a/src/python/test/test_rips_complex.py b/src/python/test/test_rips_complex.py
index b02a68e1..b86e7498 100755
--- a/src/python/test/test_rips_complex.py
+++ b/src/python/test/test_rips_complex.py
@@ -32,7 +32,7 @@ def test_rips_from_points():
     assert simplex_tree.num_simplices() == 10
     assert simplex_tree.num_vertices() == 4
 
-    assert simplex_tree.get_filtration() == [
+    assert list(simplex_tree.get_filtration()) == [
         ([0], 0.0),
         ([1], 0.0),
         ([2], 0.0),
@@ -44,6 +44,7 @@ def test_rips_from_points():
         ([1, 2], 1.4142135623730951),
         ([0, 3], 1.4142135623730951),
     ]
+
     assert simplex_tree.get_star([0]) == [
         ([0], 0.0),
         ([0, 1], 1.0),
@@ -95,7 +96,7 @@ def test_rips_from_distance_matrix():
     assert simplex_tree.num_simplices() == 10
     assert simplex_tree.num_vertices() == 4
 
-    assert simplex_tree.get_filtration() == [
+    assert list(simplex_tree.get_filtration()) == [
         ([0], 0.0),
         ([1], 0.0),
         ([2], 0.0),
@@ -107,6 +108,7 @@ def test_rips_from_distance_matrix():
         ([1, 2], 1.4142135623730951),
         ([0, 3], 1.4142135623730951),
     ]
+
     assert simplex_tree.get_star([0]) == [
         ([0], 0.0),
         ([0, 1], 1.0),
diff --git a/src/python/test/test_simplex_tree.py b/src/python/test/test_simplex_tree.py
index 1822c43b..f7848379 100755
--- a/src/python/test/test_simplex_tree.py
+++ b/src/python/test/test_simplex_tree.py
@@ -55,7 +55,7 @@ def test_insertion():
     assert st.filtration([1]) == 0.0
 
     # skeleton test
-    assert st.get_skeleton(2) == [
+    assert list(st.get_skeleton(2)) == [
         ([0, 1, 2], 4.0),
         ([0, 1], 0.0),
         ([0, 2], 4.0),
@@ -64,7 +64,7 @@ def test_insertion():
         ([1], 0.0),
         ([2], 4.0),
     ]
-    assert st.get_skeleton(1) == [
+    assert list(st.get_skeleton(1)) == [
         ([0, 1], 0.0),
         ([0, 2], 4.0),
         ([0], 0.0),
@@ -72,12 +72,12 @@ def test_insertion():
         ([1], 0.0),
         ([2], 4.0),
     ]
-    assert st.get_skeleton(0) == [([0], 0.0), ([1], 0.0), ([2], 4.0)]
+    assert list(st.get_skeleton(0)) == [([0], 0.0), ([1], 0.0), ([2], 4.0)]
 
     # remove_maximal_simplex test
     assert st.get_cofaces([0, 1, 2], 1) == []
     st.remove_maximal_simplex([0, 1, 2])
-    assert st.get_skeleton(2) == [
+    assert list(st.get_skeleton(2)) == [
         ([0, 1], 0.0),
         ([0, 2], 4.0),
         ([0], 0.0),
@@ -126,7 +126,8 @@ def test_expansion():
 
     assert st.num_vertices() == 7
     assert st.num_simplices() == 17
-    assert st.get_filtration() == [
+
+    assert list(st.get_filtration()) == [
         ([2], 0.1),
         ([3], 0.1),
         ([2, 3], 0.1),
@@ -151,7 +152,7 @@ def test_expansion():
     assert st.num_simplices() == 22
     st.initialize_filtration()
 
-    assert st.get_filtration() == [
+    assert list(st.get_filtration()) == [
         ([2], 0.1),
         ([3], 0.1),
         ([2, 3], 0.1),
@@ -248,3 +249,15 @@ def test_make_filtration_non_decreasing():
     assert st.filtration([3, 4, 5]) == 2.0
     assert st.filtration([3, 4]) == 2.0
     assert st.filtration([4, 5]) == 2.0
+
+def test_simplices_iterator():
+    st = SimplexTree()
+    
+    assert st.insert([0, 1, 2], filtration=4.0) == True
+    assert st.insert([2, 3, 4], filtration=2.0) == True
+
+    for simplex in st.get_simplices():
+        print("simplex is: ", simplex[0])
+        assert st.find(simplex[0]) == True
+        print("filtration is: ", simplex[1])
+        assert st.filtration(simplex[0]) == simplex[1]
diff --git a/src/python/test/test_tangential_complex.py b/src/python/test/test_tangential_complex.py
index e650e99c..8668a2e0 100755
--- a/src/python/test/test_tangential_complex.py
+++ b/src/python/test/test_tangential_complex.py
@@ -37,7 +37,7 @@ def test_tangential():
     assert st.num_simplices() == 6
     assert st.num_vertices() == 4
 
-    assert st.get_filtration() == [
+    assert list(st.get_filtration()) == [
         ([0], 0.0),
         ([1], 0.0),
         ([2], 0.0),
@@ -45,6 +45,7 @@ def test_tangential():
         ([3], 0.0),
         ([1, 3], 0.0),
     ]
+
     assert st.get_cofaces([0], 1) == [([0, 2], 0.0)]
 
     assert point_list[0] == tc.get_point(0)
diff --git a/src/python/test/test_time_delay.py b/src/python/test/test_time_delay.py
new file mode 100755
index 00000000..1ead9bca
--- /dev/null
+++ b/src/python/test/test_time_delay.py
@@ -0,0 +1,43 @@
+from gudhi.point_cloud.timedelay import TimeDelayEmbedding
+import numpy as np
+
+
+def test_normal():
+    # Sample array
+    ts = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
+    # Normal case.
+    prep = TimeDelayEmbedding()
+    pointclouds = prep(ts)
+    assert (pointclouds[0] == np.array([1, 2, 3])).all()
+    assert (pointclouds[1] == np.array([2, 3, 4])).all()
+    assert (pointclouds[2] == np.array([3, 4, 5])).all()
+    assert (pointclouds[3] == np.array([4, 5, 6])).all()
+    assert (pointclouds[4] == np.array([5, 6, 7])).all()
+    assert (pointclouds[5] == np.array([6, 7, 8])).all()
+    assert (pointclouds[6] == np.array([7, 8, 9])).all()
+    assert (pointclouds[7] == np.array([8, 9, 10])).all()
+    # Delay = 3
+    prep = TimeDelayEmbedding(delay=3)
+    pointclouds = prep(ts)
+    assert (pointclouds[0] == np.array([1, 4, 7])).all()
+    assert (pointclouds[1] == np.array([2, 5, 8])).all()
+    assert (pointclouds[2] == np.array([3, 6, 9])).all()
+    assert (pointclouds[3] == np.array([4, 7, 10])).all()
+    # Skip = 3
+    prep = TimeDelayEmbedding(skip=3)
+    pointclouds = prep(ts)
+    assert (pointclouds[0] == np.array([1, 2, 3])).all()
+    assert (pointclouds[1] == np.array([4, 5, 6])).all()
+    assert (pointclouds[2] == np.array([7, 8, 9])).all()
+    # Delay = 2 / Skip = 2
+    prep = TimeDelayEmbedding(delay=2, skip=2)
+    pointclouds = prep(ts)
+    assert (pointclouds[0] == np.array([1, 3, 5])).all()
+    assert (pointclouds[1] == np.array([3, 5, 7])).all()
+    assert (pointclouds[2] == np.array([5, 7, 9])).all()
+
+    # Vector series
+    ts = np.arange(0, 10).reshape(-1, 2)
+    prep = TimeDelayEmbedding(dim=4)
+    prep.fit([ts])
+    assert (prep.transform([ts])[0] == [[0, 1, 2, 3], [2, 3, 4, 5], [4, 5, 6, 7], [6, 7, 8, 9]]).all()
diff --git a/src/python/test/test_wasserstein_distance.py b/src/python/test/test_wasserstein_distance.py
index 6a6b217b..0d70e11a 100755
--- a/src/python/test/test_wasserstein_distance.py
+++ b/src/python/test/test_wasserstein_distance.py
@@ -17,7 +17,7 @@ __author__ = "Theo Lacombe"
 __copyright__ = "Copyright (C) 2019 Inria"
 __license__ = "MIT"
 
-def _basic_wasserstein(wasserstein_distance, delta, test_infinity=True):
+def _basic_wasserstein(wasserstein_distance, delta, test_infinity=True, test_matching=True):
     diag1 = np.array([[2.7, 3.7], [9.6, 14.0], [34.2, 34.974]])
     diag2 = np.array([[2.8, 4.45], [9.5, 14.1]])
     diag3 = np.array([[0, 2], [4, 6]])
@@ -51,14 +51,27 @@ def _basic_wasserstein(wasserstein_distance, delta, test_infinity=True):
     assert wasserstein_distance(diag3, diag4, internal_p=1.,     order=2.) == approx(np.sqrt(5))
     assert wasserstein_distance(diag3, diag4, internal_p=4.5,    order=2.) == approx(np.sqrt(5))
 
-    if(not test_infinity):
-        return
+    if test_infinity:
+        diag5 = np.array([[0, 3], [4, np.inf]])
+        diag6 = np.array([[7, 8], [4, 6], [3, np.inf]])
 
-    diag5 = np.array([[0, 3], [4, np.inf]])
-    diag6 = np.array([[7, 8], [4, 6], [3, np.inf]])
+        assert wasserstein_distance(diag4, diag5) == np.inf
+        assert wasserstein_distance(diag5, diag6, order=1, internal_p=np.inf) == approx(4.)
+
+
+    if test_matching:
+        match = wasserstein_distance(emptydiag, emptydiag, matching=True, internal_p=1., order=2)[1]
+        assert np.array_equal(match, [])
+        match = wasserstein_distance(emptydiag, emptydiag, matching=True, internal_p=np.inf, order=2.24)[1]
+        assert np.array_equal(match, [])
+        match = wasserstein_distance(emptydiag, diag2, matching=True, internal_p=np.inf, order=2.)[1]
+        assert np.array_equal(match , [[-1, 0], [-1, 1]])
+        match = wasserstein_distance(diag2, emptydiag, matching=True, internal_p=np.inf, order=2.24)[1]
+        assert np.array_equal(match , [[0, -1], [1, -1]])
+        match = wasserstein_distance(diag1, diag2, matching=True, internal_p=2., order=2.)[1]
+        assert np.array_equal(match, [[0, 0], [1, 1], [2, -1]])
+        
 
-    assert wasserstein_distance(diag4, diag5) == np.inf
-    assert wasserstein_distance(diag5, diag6, order=1, internal_p=np.inf) == approx(4.)
 
 def hera_wrap(delta):
     def fun(*kargs,**kwargs):
@@ -66,8 +79,8 @@ def hera_wrap(delta):
     return fun
 
 def test_wasserstein_distance_pot():
-    _basic_wasserstein(pot, 1e-15, test_infinity=False)
+    _basic_wasserstein(pot, 1e-15, test_infinity=False, test_matching=True)
 
 def test_wasserstein_distance_hera():
-    _basic_wasserstein(hera_wrap(1e-12), 1e-12)
-    _basic_wasserstein(hera_wrap(.1), .1)
+    _basic_wasserstein(hera_wrap(1e-12), 1e-12, test_matching=False)
+    _basic_wasserstein(hera_wrap(.1), .1, test_matching=False)