diff options
Diffstat (limited to 'src/python')
37 files changed, 668 insertions, 212 deletions
diff --git a/src/python/CMakeLists.txt b/src/python/CMakeLists.txt index 090a7446..f00966a5 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}") @@ -127,16 +128,6 @@ if(PYTHONINTERP_FOUND) endif () if(CGAL_FOUND) - can_cgal_use_cxx11_thread_local() - if (NOT CGAL_CAN_USE_CXX11_THREAD_LOCAL_RESULT) - if(CMAKE_BUILD_TYPE MATCHES Debug) - add_GUDHI_PYTHON_lib("${Boost_THREAD_LIBRARY_DEBUG}") - else() - add_GUDHI_PYTHON_lib("${Boost_THREAD_LIBRARY_RELEASE}") - endif() - message("** Add Boost ${Boost_LIBRARY_DIRS}") - set(GUDHI_PYTHON_LIBRARY_DIRS "${GUDHI_PYTHON_LIBRARY_DIRS}'${Boost_LIBRARY_DIRS}', ") - endif() # Add CGAL compilation args if(CGAL_HEADER_ONLY) add_gudhi_debug_info("CGAL header only version ${CGAL_VERSION}") @@ -181,7 +172,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 +217,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 +396,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_diagram_persistence_from_off_file_example.py b/src/python/example/alpha_complex_diagram_persistence_from_off_file_example.py index 4079a469..727af4fa 100755 --- a/src/python/example/alpha_complex_diagram_persistence_from_off_file_example.py +++ b/src/python/example/alpha_complex_diagram_persistence_from_off_file_example.py @@ -1,11 +1,15 @@ #!/usr/bin/env python import argparse +import errno +import os import matplotlib.pyplot as plot import gudhi -""" 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. +""" This file is part of the Gudhi Library - https://gudhi.inria.fr/ - + which is released under MIT. + See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full + license details. Author(s): Vincent Rouvreau Copyright (C) 2016 Inria @@ -41,7 +45,7 @@ args = parser.parse_args() with open(args.file, "r") as f: first_line = f.readline() if (first_line == "OFF\n") or (first_line == "nOFF\n"): - print("#####################################################################") + print("##############################################################") print("AlphaComplex creation from points read in a OFF file") message = "AlphaComplex with max_edge_length=" + repr(args.max_alpha_square) @@ -64,6 +68,7 @@ with open(args.file, "r") as f: gudhi.plot_persistence_diagram(diag, band=args.band) plot.show() else: - print(args.file, "is not a valid OFF file") + raise FileNotFoundError(errno.ENOENT, os.strerror(errno.ENOENT), + args.file) f.close() 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/alpha_rips_persistence_bottleneck_distance.py b/src/python/example/alpha_rips_persistence_bottleneck_distance.py index d5c33ec8..f156826d 100755 --- a/src/python/example/alpha_rips_persistence_bottleneck_distance.py +++ b/src/python/example/alpha_rips_persistence_bottleneck_distance.py @@ -3,9 +3,13 @@ import gudhi import argparse import math +import errno +import os -""" 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. +""" This file is part of the Gudhi Library - https://gudhi.inria.fr/ - + which is released under MIT. + See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full + license details. Author(s): Vincent Rouvreau Copyright (C) 2016 Inria @@ -36,7 +40,7 @@ with open(args.file, "r") as f: first_line = f.readline() if (first_line == "OFF\n") or (first_line == "nOFF\n"): point_cloud = gudhi.read_points_from_off_file(off_file=args.file) - print("#####################################################################") + print("##############################################################") print("RipsComplex creation from points read in a OFF file") message = "RipsComplex with max_edge_length=" + repr(args.threshold) @@ -46,14 +50,15 @@ with open(args.file, "r") as f: points=point_cloud, max_edge_length=args.threshold ) - rips_stree = rips_complex.create_simplex_tree(max_dimension=args.max_dimension) + rips_stree = rips_complex.create_simplex_tree( + max_dimension=args.max_dimension) message = "Number of simplices=" + repr(rips_stree.num_simplices()) print(message) rips_diag = rips_stree.persistence() - print("#####################################################################") + print("##############################################################") print("AlphaComplex creation from points read in a OFF file") message = "AlphaComplex with max_edge_length=" + repr(args.threshold) @@ -93,13 +98,13 @@ with open(args.file, "r") as f: print(message) max_b_distance = max(bottleneck_distance, max_b_distance) - print( - "================================================================================" - ) + print("==============================================================") message = "Bottleneck distance is " + repr(max_b_distance) print(message) else: - print(args.file, "is not a valid OFF file") + raise FileNotFoundError(errno.ENOENT, os.strerror(errno.ENOENT), + args.file) + f.close() diff --git a/src/python/example/euclidean_strong_witness_complex_diagram_persistence_from_off_file_example.py b/src/python/example/euclidean_strong_witness_complex_diagram_persistence_from_off_file_example.py index 4903667e..e1e572df 100755 --- a/src/python/example/euclidean_strong_witness_complex_diagram_persistence_from_off_file_example.py +++ b/src/python/example/euclidean_strong_witness_complex_diagram_persistence_from_off_file_example.py @@ -1,11 +1,15 @@ #!/usr/bin/env python import argparse +import errno +import os import matplotlib.pyplot as plot import gudhi -""" 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. +""" This file is part of the Gudhi Library - https://gudhi.inria.fr/ - + which is released under MIT. + See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full + license details. Author(s): Vincent Rouvreau Copyright (C) 2016 Inria @@ -44,8 +48,9 @@ args = parser.parse_args() with open(args.file, "r") as f: first_line = f.readline() if (first_line == "OFF\n") or (first_line == "nOFF\n"): - print("#####################################################################") - print("EuclideanStrongWitnessComplex creation from points read in a OFF file") + print("##############################################################") + print("EuclideanStrongWitnessComplex creation from points read "\ + "in a OFF file") witnesses = gudhi.read_points_from_off_file(off_file=args.file) landmarks = gudhi.pick_n_random_points( @@ -64,7 +69,8 @@ with open(args.file, "r") as f: witnesses=witnesses, landmarks=landmarks ) simplex_tree = witness_complex.create_simplex_tree( - max_alpha_square=args.max_alpha_square, limit_dimension=args.limit_dimension + max_alpha_square=args.max_alpha_square, + limit_dimension=args.limit_dimension ) message = "Number of simplices=" + repr(simplex_tree.num_simplices()) @@ -79,6 +85,7 @@ with open(args.file, "r") as f: gudhi.plot_persistence_diagram(diag, band=args.band) plot.show() else: - print(args.file, "is not a valid OFF file") + raise FileNotFoundError(errno.ENOENT, os.strerror(errno.ENOENT), + args.file) f.close() diff --git a/src/python/example/euclidean_witness_complex_diagram_persistence_from_off_file_example.py b/src/python/example/euclidean_witness_complex_diagram_persistence_from_off_file_example.py index 339a8577..58cb2bb5 100755 --- a/src/python/example/euclidean_witness_complex_diagram_persistence_from_off_file_example.py +++ b/src/python/example/euclidean_witness_complex_diagram_persistence_from_off_file_example.py @@ -1,11 +1,15 @@ #!/usr/bin/env python import argparse +import errno +import os import matplotlib.pyplot as plot import gudhi -""" 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. +""" This file is part of the Gudhi Library - https://gudhi.inria.fr/ - + which is released under MIT. + See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full + license details. Author(s): Vincent Rouvreau Copyright (C) 2016 Inria @@ -78,6 +82,7 @@ with open(args.file, "r") as f: gudhi.plot_persistence_diagram(diag, band=args.band) plot.show() else: - print(args.file, "is not a valid OFF file") + raise FileNotFoundError(errno.ENOENT, os.strerror(errno.ENOENT), + args.file) f.close() diff --git a/src/python/example/periodic_cubical_complex_barcode_persistence_from_perseus_file_example.py b/src/python/example/periodic_cubical_complex_barcode_persistence_from_perseus_file_example.py index c692e66f..499171df 100755 --- a/src/python/example/periodic_cubical_complex_barcode_persistence_from_perseus_file_example.py +++ b/src/python/example/periodic_cubical_complex_barcode_persistence_from_perseus_file_example.py @@ -2,10 +2,14 @@ import argparse import matplotlib.pyplot as plot +import errno +import os import gudhi -""" 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. +""" This file is part of the Gudhi Library - https://gudhi.inria.fr/ - + which is released under MIT. + See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full + license details. Author(s): Vincent Rouvreau Copyright (C) 2016 Inria @@ -57,9 +61,10 @@ parser.add_argument( args = parser.parse_args() if is_file_perseus(args.file): - print("#####################################################################") + print("##################################################################") print("PeriodicCubicalComplex creation") - periodic_cubical_complex = gudhi.PeriodicCubicalComplex(perseus_file=args.file) + periodic_cubical_complex = gudhi.PeriodicCubicalComplex( + perseus_file=args.file) print("persistence(homology_coeff_field=3, min_persistence=0)=") diag = periodic_cubical_complex.persistence( @@ -73,4 +78,5 @@ if is_file_perseus(args.file): gudhi.plot_persistence_barcode(diag) plot.show() else: - print(args.file, "is not a valid perseus style file") + raise FileNotFoundError(errno.ENOENT, os.strerror(errno.ENOENT), + args.file) diff --git a/src/python/example/rips_complex_diagram_persistence_from_off_file_example.py b/src/python/example/rips_complex_diagram_persistence_from_off_file_example.py index c757aca7..6f992508 100755 --- a/src/python/example/rips_complex_diagram_persistence_from_off_file_example.py +++ b/src/python/example/rips_complex_diagram_persistence_from_off_file_example.py @@ -1,11 +1,15 @@ #!/usr/bin/env python import argparse +import errno +import os import matplotlib.pyplot as plot import gudhi -""" 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. +""" This file is part of the Gudhi Library - https://gudhi.inria.fr/ - + which is released under MIT. + See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full + license details. Author(s): Vincent Rouvreau Copyright (C) 2016 Inria @@ -42,10 +46,11 @@ args = parser.parse_args() with open(args.file, "r") as f: first_line = f.readline() if (first_line == "OFF\n") or (first_line == "nOFF\n"): - print("#####################################################################") + print("##############################################################") print("RipsComplex creation from points read in a OFF file") - message = "RipsComplex with max_edge_length=" + repr(args.max_edge_length) + message = "RipsComplex with max_edge_length=" + \ + repr(args.max_edge_length) print(message) point_cloud = gudhi.read_points_from_off_file(off_file=args.file) @@ -68,6 +73,7 @@ with open(args.file, "r") as f: gudhi.plot_persistence_diagram(diag, band=args.band) plot.show() else: - print(args.file, "is not a valid OFF file") + raise FileNotFoundError(errno.ENOENT, os.strerror(errno.ENOENT), + args.file) f.close() 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/example/tangential_complex_plain_homology_from_off_file_example.py b/src/python/example/tangential_complex_plain_homology_from_off_file_example.py index f0df2189..85bade4a 100755 --- a/src/python/example/tangential_complex_plain_homology_from_off_file_example.py +++ b/src/python/example/tangential_complex_plain_homology_from_off_file_example.py @@ -1,11 +1,15 @@ #!/usr/bin/env python import argparse +import errno +import os import matplotlib.pyplot as plot import gudhi -""" 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. +""" This file is part of the Gudhi Library - https://gudhi.inria.fr/ - + which is released under MIT. + See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full + license details. Author(s): Vincent Rouvreau Copyright (C) 2016 Inria @@ -19,7 +23,7 @@ __copyright__ = "Copyright (C) 2016 Inria" __license__ = "MIT" parser = argparse.ArgumentParser( - description="TangentialComplex creation from " "points read in a OFF file.", + description="TangentialComplex creation from points read in a OFF file.", epilog="Example: " "example/tangential_complex_plain_homology_from_off_file_example.py " "-f ../data/points/tore3D_300.off -i 3" @@ -41,10 +45,11 @@ args = parser.parse_args() with open(args.file, "r") as f: first_line = f.readline() if (first_line == "OFF\n") or (first_line == "nOFF\n"): - print("#####################################################################") + print("##############################################################") print("TangentialComplex creation from points read in a OFF file") - tc = gudhi.TangentialComplex(intrisic_dim=args.intrisic_dim, off_file=args.file) + tc = gudhi.TangentialComplex(intrisic_dim=args.intrisic_dim, + off_file=args.file) tc.compute_tangential_complex() st = tc.create_simplex_tree() @@ -60,6 +65,7 @@ with open(args.file, "r") as f: gudhi.plot_persistence_diagram(diag, band=args.band) plot.show() else: - print(args.file, "is not a valid OFF file") + raise FileNotFoundError(errno.ENOENT, os.strerror(errno.ENOENT), + args.file) f.close() diff --git a/src/python/gudhi/alpha_complex.pyx b/src/python/gudhi/alpha_complex.pyx index fff3e920..e04dc652 100644 --- a/src/python/gudhi/alpha_complex.pyx +++ b/src/python/gudhi/alpha_complex.pyx @@ -1,5 +1,7 @@ -# 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. +# This file is part of the Gudhi Library - https://gudhi.inria.fr/ - +# which is released under MIT. +# See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full +# license details. # Author(s): Vincent Rouvreau # # Copyright (C) 2016 Inria @@ -7,6 +9,7 @@ # Modification(s): # - YYYY/MM Author: Description of the modification +from __future__ import print_function from cython cimport numeric from libcpp.vector cimport vector from libcpp.utility cimport pair @@ -69,7 +72,8 @@ cdef class AlphaComplex: def __cinit__(self, points = None, off_file = ''): if off_file: if os.path.isfile(off_file): - self.thisptr = new Alpha_complex_interface(off_file.encode('utf-8'), True) + self.thisptr = new Alpha_complex_interface( + off_file.encode('utf-8'), True) else: print("file " + off_file + " not found.") else: diff --git a/src/python/gudhi/cubical_complex.pyx b/src/python/gudhi/cubical_complex.pyx index cbeda014..d5ad1266 100644 --- a/src/python/gudhi/cubical_complex.pyx +++ b/src/python/gudhi/cubical_complex.pyx @@ -1,5 +1,7 @@ -# 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. +# This file is part of the Gudhi Library - https://gudhi.inria.fr/ - +# which is released under MIT. +# See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full +# license details. # Author(s): Vincent Rouvreau # # Copyright (C) 2016 Inria @@ -7,12 +9,15 @@ # Modification(s): # - YYYY/MM Author: Description of the modification +from __future__ import print_function from cython cimport numeric from libcpp.vector cimport vector from libcpp.utility cimport pair from libcpp.string cimport string from libcpp cimport bool +import errno import os +import sys import numpy as np @@ -87,10 +92,12 @@ cdef class CubicalComplex: if os.path.isfile(perseus_file): self.thisptr = new Bitmap_cubical_complex_base_interface(perseus_file.encode('utf-8')) else: - print("file " + perseus_file + " not found.") + raise FileNotFoundError(errno.ENOENT, os.strerror(errno.ENOENT), + perseus_file) else: print("CubicalComplex can be constructed from dimensions and " - "top_dimensional_cells or from a Perseus-style file name.") + "top_dimensional_cells or from a Perseus-style file name.", + file=sys.stderr) def __dealloc__(self): if self.thisptr != NULL: @@ -199,5 +206,5 @@ cdef class CubicalComplex: intervals_result = self.pcohptr.intervals_in_dimension(dimension) else: print("intervals_in_dim function requires persistence function" - " to be launched first.") + " to be launched first.", file=sys.stderr) return np.array(intervals_result) diff --git a/src/python/gudhi/nerve_gic.pyx b/src/python/gudhi/nerve_gic.pyx index 45cc8eba..9c89b239 100644 --- a/src/python/gudhi/nerve_gic.pyx +++ b/src/python/gudhi/nerve_gic.pyx @@ -1,5 +1,7 @@ -# 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. +# This file is part of the Gudhi Library - https://gudhi.inria.fr/ - +# which is released under MIT. +# See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full +# license details. # Author(s): Vincent Rouvreau # # Copyright (C) 2018 Inria @@ -7,11 +9,13 @@ # Modification(s): # - YYYY/MM Author: Description of the modification +from __future__ import print_function from cython cimport numeric from libcpp.vector cimport vector from libcpp.utility cimport pair from libcpp.string cimport string from libcpp cimport bool +import errno import os from libc.stdint cimport intptr_t @@ -96,7 +100,8 @@ cdef class CoverComplex: return self.thisptr != NULL def set_point_cloud_from_range(self, cloud): - """ Reads and stores the input point cloud from a vector stored in memory. + """ Reads and stores the input point cloud from a vector stored in + memory. :param cloud: Input vector containing the point cloud. :type cloud: vector[vector[double]] @@ -104,7 +109,8 @@ cdef class CoverComplex: return self.thisptr.set_point_cloud_from_range(cloud) def set_distances_from_range(self, distance_matrix): - """ Reads and stores the input distance matrix from a vector stored in memory. + """ Reads and stores the input distance matrix from a vector stored in + memory. :param distance_matrix: Input vector containing the distance matrix. :type distance_matrix: vector[vector[double]] @@ -163,7 +169,8 @@ cdef class CoverComplex: """ stree = SimplexTree() cdef intptr_t stree_int_ptr=stree.thisptr - self.thisptr.create_simplex_tree(<Simplex_tree_interface_full_featured*>stree_int_ptr) + self.thisptr.create_simplex_tree( + <Simplex_tree_interface_full_featured*>stree_int_ptr) return stree def find_simplices(self): @@ -182,8 +189,8 @@ cdef class CoverComplex: if os.path.isfile(off_file): return self.thisptr.read_point_cloud(off_file.encode('utf-8')) else: - print("file " + off_file + " not found.") - return False + raise FileNotFoundError(errno.ENOENT, os.strerror(errno.ENOENT), + off_file) def set_automatic_resolution(self): """Computes the optimal length of intervals (i.e. the smallest interval @@ -214,7 +221,8 @@ cdef class CoverComplex: if os.path.isfile(color_file_name): self.thisptr.set_color_from_file(color_file_name.encode('utf-8')) else: - print("file " + color_file_name + " not found.") + raise FileNotFoundError(errno.ENOENT, os.strerror(errno.ENOENT), + color_file_name) def set_color_from_range(self, color): """Computes the function used to color the nodes of the simplicial @@ -235,7 +243,8 @@ cdef class CoverComplex: if os.path.isfile(cover_file_name): self.thisptr.set_cover_from_file(cover_file_name.encode('utf-8')) else: - print("file " + cover_file_name + " not found.") + raise FileNotFoundError(errno.ENOENT, os.strerror(errno.ENOENT), + cover_file_name) def set_cover_from_function(self): """Creates a cover C from the preimages of the function f. @@ -268,7 +277,8 @@ cdef class CoverComplex: if os.path.isfile(func_file_name): self.thisptr.set_function_from_file(func_file_name.encode('utf-8')) else: - print("file " + func_file_name + " not found.") + raise FileNotFoundError(errno.ENOENT, os.strerror(errno.ENOENT), + func_file_name) def set_function_from_range(self, function): """Creates the function f from a vector stored in memory. @@ -302,14 +312,15 @@ cdef class CoverComplex: """Creates a graph G from a file containing the edges. :param graph_file_name: Name of the input graph file. The graph file - contains one edge per line, each edge being represented by the IDs of - its two nodes. + contains one edge per line, each edge being represented by the IDs + of its two nodes. :type graph_file_name: string """ if os.path.isfile(graph_file_name): self.thisptr.set_graph_from_file(graph_file_name.encode('utf-8')) else: - print("file " + graph_file_name + " not found.") + raise FileNotFoundError(errno.ENOENT, os.strerror(errno.ENOENT), + graph_file_name) def set_graph_from_OFF(self): """Creates a graph G from the triangulation given by the input OFF diff --git a/src/python/gudhi/off_reader.pyx b/src/python/gudhi/off_reader.pyx index 7e6d9d80..a3200704 100644 --- a/src/python/gudhi/off_reader.pyx +++ b/src/python/gudhi/off_reader.pyx @@ -1,5 +1,7 @@ -# 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. +# This file is part of the Gudhi Library - https://gudhi.inria.fr/ - +# which is released under MIT. +# See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full +# license details. # Author(s): Vincent Rouvreau # # Copyright (C) 2016 Inria @@ -7,9 +9,11 @@ # Modification(s): # - YYYY/MM Author: Description of the modification +from __future__ import print_function from cython cimport numeric from libcpp.vector cimport vector from libcpp.string cimport string +import errno import os __author__ = "Vincent Rouvreau" @@ -32,6 +36,6 @@ def read_points_from_off_file(off_file=''): if os.path.isfile(off_file): return read_points_from_OFF_file(off_file.encode('utf-8')) else: - print("file " + off_file + " not found.") - return [] + raise FileNotFoundError(errno.ENOENT, os.strerror(errno.ENOENT), + off_file) diff --git a/src/python/gudhi/periodic_cubical_complex.pyx b/src/python/gudhi/periodic_cubical_complex.pyx index 37f76201..fd08b976 100644 --- a/src/python/gudhi/periodic_cubical_complex.pyx +++ b/src/python/gudhi/periodic_cubical_complex.pyx @@ -7,11 +7,13 @@ # Modification(s): # - YYYY/MM Author: Description of the modification +from __future__ import print_function from cython cimport numeric from libcpp.vector cimport vector from libcpp.utility cimport pair from libcpp.string cimport string from libcpp cimport bool +import sys import os import numpy as np @@ -95,12 +97,12 @@ cdef class PeriodicCubicalComplex: if os.path.isfile(perseus_file): self.thisptr = new Periodic_cubical_complex_base_interface(perseus_file.encode('utf-8')) else: - print("file " + perseus_file + " not found.") + print("file " + perseus_file + " not found.", file=sys.stderr) else: print("CubicalComplex can be constructed from dimensions, " "top_dimensional_cells and periodic_dimensions, or from " "top_dimensional_cells and periodic_dimensions or from " - "a Perseus-style file name.") + "a Perseus-style file name.", file=sys.stderr) def __dealloc__(self): if self.thisptr != NULL: @@ -209,5 +211,5 @@ cdef class PeriodicCubicalComplex: intervals_result = self.pcohptr.intervals_in_dimension(dimension) else: print("intervals_in_dim function requires persistence function" - " to be launched first.") + " to be launched first.", file=sys.stderr) return np.array(intervals_result) 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 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_cover_complex.py b/src/python/test/test_cover_complex.py index 32bc5a26..260f6a5c 100755 --- a/src/python/test/test_cover_complex.py +++ b/src/python/test/test_cover_complex.py @@ -9,6 +9,7 @@ """ from gudhi import CoverComplex +import pytest __author__ = "Vincent Rouvreau" __copyright__ = "Copyright (C) 2018 Inria" @@ -24,7 +25,8 @@ def test_empty_constructor(): def test_non_existing_file_read(): # Try to open a non existing file cover = CoverComplex() - assert cover.read_point_cloud("pouetpouettralala.toubiloubabdou") == False + with pytest.raises(FileNotFoundError): + cover.read_point_cloud("pouetpouettralala.toubiloubabdou") def test_files_creation(): diff --git a/src/python/test/test_cubical_complex.py b/src/python/test/test_cubical_complex.py index 8c1b2600..fce4875c 100755 --- a/src/python/test/test_cubical_complex.py +++ b/src/python/test/test_cubical_complex.py @@ -10,6 +10,7 @@ from gudhi import CubicalComplex, PeriodicCubicalComplex import numpy as np +import pytest __author__ = "Vincent Rouvreau" __copyright__ = "Copyright (C) 2016 Inria" @@ -25,9 +26,8 @@ def test_empty_constructor(): def test_non_existing_perseus_file_constructor(): # Try to open a non existing file - cub = CubicalComplex(perseus_file="pouetpouettralala.toubiloubabdou") - assert cub.__is_defined() == False - assert cub.__is_persistence_defined() == False + with pytest.raises(FileNotFoundError): + cub = CubicalComplex(perseus_file="pouetpouettralala.toubiloubabdou") def test_dimension_or_perseus_file_constructor(): 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_subsampling.py b/src/python/test/test_subsampling.py index fe0985fa..31f64e32 100755 --- a/src/python/test/test_subsampling.py +++ b/src/python/test/test_subsampling.py @@ -120,7 +120,6 @@ def test_simple_pick_n_random_points(): # Go furter than point set on purpose for iter in range(1, 10): sub_set = gudhi.pick_n_random_points(points=point_set, nb_points=iter) - print(5) for sub in sub_set: found = False for point in point_set: 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) |