Merge remote-tracking branch 'origin/master' into gen2

22 files changed, 528 insertions, 131 deletions

diff --git a/src/python/CMakeLists.txt b/src/python/CMakeLists.txt
index 090a7446..22af3ec9 100644
--- a/src/python/CMakeLists.txt
+++ b/src/python/CMakeLists.txt
@@ -56,6 +56,7 @@ if(PYTHONINTERP_FOUND)
     # Modules that should not be auto-imported in __init__.py
     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}'point_cloud', ")
 
     add_gudhi_debug_info("Python version ${PYTHON_VERSION_STRING}")
     add_gudhi_debug_info("Cython version ${CYTHON_VERSION}")
@@ -181,7 +182,7 @@ if(PYTHONINTERP_FOUND)
         set(GUDHI_PYTHON_LIBRARY_DIRS "${GUDHI_PYTHON_LIBRARY_DIRS}'${MPFR_LIBRARIES_DIR}', ")
         message("** Add mpfr ${MPFR_LIBRARIES}")
       endif(MPFR_FOUND)
-endif(CGAL_FOUND)
+    endif(CGAL_FOUND)
 
     # Specific for Mac
     if (${CMAKE_SYSTEM_NAME} MATCHES "Darwin")
@@ -226,6 +227,7 @@ endif(CGAL_FOUND)
     file(COPY "gudhi/persistence_graphical_tools.py" DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/gudhi")
     file(COPY "gudhi/representations" DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/gudhi/")
     file(COPY "gudhi/wasserstein.py" DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/gudhi")
+    file(COPY "gudhi/point_cloud" DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/gudhi")
 
     add_custom_command(
         OUTPUT gudhi.so
@@ -404,6 +406,9 @@ endif(CGAL_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/persistence_graphical_tools_sum.inc b/src/python/doc/persistence_graphical_tools_sum.inc
index 0cdf8072..ef376802 100644
--- a/src/python/doc/persistence_graphical_tools_sum.inc
+++ b/src/python/doc/persistence_graphical_tools_sum.inc
@@ -2,11 +2,11 @@
    :widths: 30 50 20
 
    +-----------------------------------------------------------------+-----------------------------------------------------------------------+-----------------------------------------------+
-   | .. figure::                                                     | These graphical tools comes on top of persistence results and allows  | :Author: Vincent Rouvreau                     |
-   |      img/graphical_tools_representation.png                     | the user to build easily persistence barcode, diagram or density.     |                                               |
+   | .. figure::                                                     | These graphical tools comes on top of persistence results and allows  | :Author: Vincent Rouvreau, Theo Lacombe       |
+   |      img/graphical_tools_representation.png                     | the user to display easily persistence barcode, diagram or density.   |                                               |
    |                                                                 |                                                                       | :Introduced in: GUDHI 2.0.0                   |
-   |                                                                 |                                                                       |                                               |
-   |                                                                 |                                                                       | :Copyright: MIT                               |
+   |                                                                 | Note that these functions return the matplotlib axis, allowing        |                                               |
+   |                                                                 | for further modifications (title, aspect, etc.)                       | :Copyright: MIT                               |
    |                                                                 |                                                                       |                                               |
    |                                                                 |                                                                       | :Requires: matplotlib, numpy and scipy        |
    +-----------------------------------------------------------------+-----------------------------------------------------------------------+-----------------------------------------------+
diff --git a/src/python/doc/persistence_graphical_tools_user.rst b/src/python/doc/persistence_graphical_tools_user.rst
index 80002db6..91e52703 100644
--- a/src/python/doc/persistence_graphical_tools_user.rst
+++ b/src/python/doc/persistence_graphical_tools_user.rst
@@ -20,7 +20,7 @@ This function can display the persistence result as a barcode:
 .. plot::
    :include-source:
 
-    import matplotlib.pyplot as plot
+    import matplotlib.pyplot as plt
     import gudhi
 
     off_file = gudhi.__root_source_dir__ + '/data/points/tore3D_300.off'
@@ -31,7 +31,7 @@ This function can display the persistence result as a barcode:
     diag = simplex_tree.persistence(min_persistence=0.4)
 
     gudhi.plot_persistence_barcode(diag)
-    plot.show()
+    plt.show()
 
 Show persistence as a diagram
 -----------------------------
@@ -44,15 +44,31 @@ This function can display the persistence result as a diagram:
 .. plot::
    :include-source:
 
-    import matplotlib.pyplot as plot
+    import matplotlib.pyplot as plt
     import gudhi
 
     # rips_on_tore3D_1307.pers obtained from write_persistence_diagram method
     persistence_file=gudhi.__root_source_dir__ + \
         '/data/persistence_diagram/rips_on_tore3D_1307.pers'
-    gudhi.plot_persistence_diagram(persistence_file=persistence_file,
+    ax = gudhi.plot_persistence_diagram(persistence_file=persistence_file,
         legend=True)
-    plot.show()
+    # We can modify the title, aspect, etc.
+    ax.set_title("Persistence diagram of a torus")
+    ax.set_aspect("equal")  # forces to be square shaped
+    plt.show()
+
+Note that (as barcode and density) it can also take a simple `np.array` 
+of shape (N x 2) encoding a persistence diagram (in a given dimension).
+
+.. plot::
+    :include-source:
+    
+    import matplotlib.pyplot as plt
+    import gudhi
+    import numpy as np
+    d = np.array([[0, 1], [1, 2], [1, np.inf]])
+    gudhi.plot_persistence_diagram(d)
+    plt.show()
 
 Persistence density
 -------------------
@@ -65,7 +81,7 @@ If you want more information on a specific dimension, for instance:
 .. plot::
    :include-source:
 
-    import matplotlib.pyplot as plot
+    import matplotlib.pyplot as plt
     import gudhi
     # rips_on_tore3D_1307.pers obtained from write_persistence_diagram method
     persistence_file=gudhi.__root_source_dir__ + \
@@ -75,9 +91,9 @@ If you want more information on a specific dimension, for instance:
         only_this_dim=1)
     pers_diag = [(1, elt) for elt in birth_death]
     # Use subplots to display diagram and density side by side
-    fig, axes = plot.subplots(nrows=1, ncols=2, figsize=(12, 5))
+    fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(12, 5))
     gudhi.plot_persistence_diagram(persistence=pers_diag,
         axes=axes[0])
     gudhi.plot_persistence_density(persistence=pers_diag,
         dimension=1, legend=True, axes=axes[1])
-    plot.show()
+    plt.show()
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/persistence_graphical_tools.py b/src/python/gudhi/persistence_graphical_tools.py
index 246280de..cc3db467 100644
--- a/src/python/gudhi/persistence_graphical_tools.py
+++ b/src/python/gudhi/persistence_graphical_tools.py
@@ -5,6 +5,7 @@
 # Copyright (C) 2016 Inria
 #
 # Modification(s):
+#   - 2020/02 Theo Lacombe: Added more options for improved rendering and more flexibility.
 #   - YYYY/MM Author: Description of the modification
 
 from os import path
@@ -14,7 +15,7 @@ import numpy as np
 from gudhi.reader_utils import read_persistence_intervals_in_dimension
 from gudhi.reader_utils import read_persistence_intervals_grouped_by_dimension
 
-__author__ = "Vincent Rouvreau, Bertrand Michel"
+__author__ = "Vincent Rouvreau, Bertrand Michel, Theo Lacombe"
 __copyright__ = "Copyright (C) 2016 Inria"
 __license__ = "MIT"
 
@@ -43,6 +44,19 @@ def __min_birth_max_death(persistence, band=0.0):
         max_death += band
     return (min_birth, max_death)
 
+
+def _array_handler(a):
+    '''
+    :param a: if array, assumes it is a (n x 2) np.array and return a
+                persistence-compatible list (padding with 0), so that the
+                plot can be performed seamlessly.
+    '''
+    if isinstance(a[0][1], np.float64) or isinstance(a[0][1], float):
+        return [[0, x] for x in a]
+    else:
+        return a
+
+
 def plot_persistence_barcode(
     persistence=[],
     persistence_file="",
@@ -52,13 +66,16 @@ def plot_persistence_barcode(
     inf_delta=0.1,
     legend=False,
     colormap=None,
-    axes=None
+    axes=None,
+    fontsize=16,
 ):
     """This function plots the persistence bar code from persistence values list
+    , a np.array of shape (N x 2) (representing a diagram 
+    in a single homology dimension), 
     or from a :doc:`persistence file <fileformats>`.
 
-    :param persistence: Persistence intervals values list grouped by dimension.
-    :type persistence: list of tuples(dimension, tuple(birth, death)).
+    :param persistence: Persistence intervals values list. Can be grouped by dimension or not.
+    :type persistence: an array of (dimension, array of (birth, death)) or an array of (birth, death).
     :param persistence_file: A :doc:`persistence file <fileformats>` style name
         (reset persistence if both are set).
     :type persistence_file: string
@@ -81,11 +98,19 @@ def plot_persistence_barcode(
     :param axes: A matplotlib-like subplot axes. If None, the plot is drawn on
         a new set of axes.
     :type axes: `matplotlib.axes.Axes`
+    :param fontsize: Fontsize to use in axis.
+    :type fontsize: int
     :returns: (`matplotlib.axes.Axes`): The axes on which the plot was drawn.
     """
     try:
         import matplotlib.pyplot as plt
         import matplotlib.patches as mpatches
+        from matplotlib import rc
+        plt.rc('text', usetex=True)
+        plt.rc('font', family='serif')
+
+
+        persistence = _array_handler(persistence)
 
         if persistence_file != "":
             if path.isfile(persistence_file):
@@ -163,7 +188,7 @@ def plot_persistence_barcode(
                 loc="lower right",
             )
 
-        axes.set_title("Persistence barcode")
+        axes.set_title("Persistence barcode", fontsize=fontsize)
 
         # Ends plot on infinity value and starts a little bit before min_birth
         axes.axis([axis_start, infinity, 0, ind])
@@ -183,13 +208,16 @@ def plot_persistence_diagram(
     inf_delta=0.1,
     legend=False,
     colormap=None,
-    axes=None
+    axes=None,
+    fontsize=16,
+    greyblock=True
 ):
     """This function plots the persistence diagram from persistence values
-    list or from a :doc:`persistence file <fileformats>`.
+    list, a np.array of shape (N x 2) representing a diagram in a single
+    homology dimension, or from a :doc:`persistence file <fileformats>`.
 
-    :param persistence: Persistence intervals values list grouped by dimension.
-    :type persistence: list of tuples(dimension, tuple(birth, death)).
+    :param persistence: Persistence intervals values list. Can be grouped by dimension or not.
+    :type persistence: an array of (dimension, array of (birth, death)) or an array of (birth, death).
     :param persistence_file: A :doc:`persistence file <fileformats>` style name
         (reset persistence if both are set).
     :type persistence_file: string
@@ -214,11 +242,20 @@ def plot_persistence_diagram(
     :param axes: A matplotlib-like subplot axes. If None, the plot is drawn on
         a new set of axes.
     :type axes: `matplotlib.axes.Axes`
+    :param fontsize: Fontsize to use in axis.
+    :type fontsize: int
+    :param greyblock: if we want to plot a grey patch on the lower half plane for nicer rendering. Default True.
+    :type greyblock: boolean
     :returns: (`matplotlib.axes.Axes`): The axes on which the plot was drawn.
     """
     try:
         import matplotlib.pyplot as plt
         import matplotlib.patches as mpatches
+        from matplotlib import rc
+        plt.rc('text', usetex=True)
+        plt.rc('font', family='serif')
+
+        persistence = _array_handler(persistence)
 
         if persistence_file != "":
             if path.isfile(persistence_file):
@@ -256,19 +293,18 @@ def plot_persistence_diagram(
         # Replace infinity values with max_death + delta for diagram to be more
         # readable
         infinity = max_death + delta
+        axis_end = max_death + delta / 2
         axis_start = min_birth - delta
 
-        # line display of equation : birth = death
-        x = np.linspace(axis_start, infinity, 1000)
-        # infinity line and text
-        axes.plot(x, x, color="k", linewidth=1.0)
-        axes.plot(x, [infinity] * len(x), linewidth=1.0, color="k", alpha=alpha)
-        axes.text(axis_start, infinity, r"$\infty$", color="k", alpha=alpha)
         # bootstrap band
         if band > 0.0:
+            x = np.linspace(axis_start, infinity, 1000)
             axes.fill_between(x, x, x + band, alpha=alpha, facecolor="red")
-
+        # lower diag patch
+        if greyblock:
+            axes.add_patch(mpatches.Polygon([[axis_start, axis_start], [axis_end, axis_start], [axis_end, axis_end]], fill=True, color='lightgrey'))
         # Draw points in loop
+        pts_at_infty = False  # Records presence of pts at infty
         for interval in reversed(persistence):
             if float(interval[1][1]) != float("inf"):
                 # Finite death case
@@ -279,10 +315,23 @@ def plot_persistence_diagram(
                     color=colormap[interval[0]],
                 )
             else:
+                pts_at_infty = True
                 # Infinite death case for diagram to be nicer
                 axes.scatter(
                     interval[1][0], infinity, alpha=alpha, color=colormap[interval[0]]
                 )
+        if pts_at_infty:
+            # infinity line and text
+            axes.plot([axis_start, axis_end], [axis_start, axis_end], linewidth=1.0, color="k")
+            axes.plot([axis_start, axis_end], [infinity, infinity], linewidth=1.0, color="k", alpha=alpha)
+            # Infinity label
+            yt = axes.get_yticks()
+            yt = yt[np.where(yt < axis_end)] # to avoid ploting ticklabel higher than infinity
+            yt = np.append(yt, infinity)
+            ytl = ["%.3f" % e for e in yt]  # to avoid float precision error
+            ytl[-1] = r'$+\infty$'
+            axes.set_yticks(yt)
+            axes.set_yticklabels(ytl)
 
         if legend:
             dimensions = list(set(item[0] for item in persistence))
@@ -293,11 +342,11 @@ def plot_persistence_diagram(
                 ]
             )
 
-        axes.set_xlabel("Birth")
-        axes.set_ylabel("Death")
+        axes.set_xlabel("Birth", fontsize=fontsize)
+        axes.set_ylabel("Death", fontsize=fontsize)
+        axes.set_title("Persistence diagram", fontsize=fontsize)
         # Ends plot on infinity value and starts a little bit before min_birth
-        axes.axis([axis_start, infinity, axis_start, infinity + delta])
-        axes.set_title("Persistence diagram")
+        axes.axis([axis_start, axis_end, axis_start, infinity + delta/2])
         return axes
 
     except ImportError:
@@ -313,16 +362,22 @@ def plot_persistence_density(
     dimension=None,
     cmap=None,
     legend=False,
-    axes=None
+    axes=None,
+    fontsize=16,
+    greyblock=False
 ):
     """This function plots the persistence density from persistence
-    values list or from a :doc:`persistence file <fileformats>`. Be
+    values list, np.array of shape (N x 2) representing a diagram
+    in a single homology dimension,
+    or from a :doc:`persistence file <fileformats>`. Be
     aware that this function does not distinguish the dimension, it is
     up to you to select the required one. This function also does not handle
     degenerate data set (scipy correlation matrix inversion can fail).
 
-    :param persistence: Persistence intervals values list grouped by dimension.
-    :type persistence: list of tuples(dimension, tuple(birth, death)).
+    :param persistence: Persistence intervals values list. 
+                        Can be grouped by dimension or not.
+    :type persistence: an array of (dimension, array of (birth, death)) 
+                        or an array of (birth, death).
     :param persistence_file: A :doc:`persistence file <fileformats>`
         style name (reset persistence if both are set).
     :type persistence_file: string
@@ -355,11 +410,22 @@ def plot_persistence_density(
     :param axes: A matplotlib-like subplot axes. If None, the plot is drawn on
         a new set of axes.
     :type axes: `matplotlib.axes.Axes`
+    :param fontsize: Fontsize to use in axis.
+    :type fontsize: int
+    :param greyblock: if we want to plot a grey patch on the lower half plane 
+                         for nicer rendering. Default False.
+    :type greyblock: boolean
     :returns: (`matplotlib.axes.Axes`): The axes on which the plot was drawn.
     """
     try:
         import matplotlib.pyplot as plt
+        import matplotlib.patches as mpatches
         from scipy.stats import kde
+        from matplotlib import rc
+        plt.rc('text', usetex=True)
+        plt.rc('font', family='serif')
+
+        persistence = _array_handler(persistence)
 
         if persistence_file != "":
             if dimension is None:
@@ -418,12 +484,16 @@ def plot_persistence_density(
         # Make the plot
         img = axes.pcolormesh(xi, yi, zi.reshape(xi.shape), cmap=cmap)
 
+        if greyblock:
+            axes.add_patch(mpatches.Polygon([[birth.min(), birth.min()], [death.max(), birth.min()], [death.max(), death.max()]], fill=True, color='lightgrey'))
+
         if legend:
             plt.colorbar(img, ax=axes)
 
-        axes.set_xlabel("Birth")
-        axes.set_ylabel("Death")
-        axes.set_title("Persistence density")
+        axes.set_xlabel("Birth", fontsize=fontsize)
+        axes.set_ylabel("Death", fontsize=fontsize)
+        axes.set_title("Persistence density", fontsize=fontsize)
+
         return axes
 
     except ImportError:
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 53e2bbc9..44789365 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 d5f642d1..faa9f9d8 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
 import numpy
 from numpy import array as np_array
@@ -208,22 +209,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
@@ -234,15 +247,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)