From a064f5698fedbe13f6c343cb0b82e0f4d72caffb Mon Sep 17 00:00:00 2001 From: ROUVREAU Vincent Date: Mon, 27 Jan 2020 17:37:31 +0100 Subject: A first naive iterator implementation with yield --- src/python/gudhi/simplex_tree.pxd | 8 +++++++- src/python/gudhi/simplex_tree.pyx | 18 +++++++++--------- src/python/include/Simplex_tree_interface.h | 28 +++++++++++++++++----------- 3 files changed, 33 insertions(+), 21 deletions(-) (limited to 'src/python') diff --git a/src/python/gudhi/simplex_tree.pxd b/src/python/gudhi/simplex_tree.pxd index 96d14079..caf3c459 100644 --- a/src/python/gudhi/simplex_tree.pxd +++ b/src/python/gudhi/simplex_tree.pxd @@ -21,6 +21,9 @@ 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::Simplex_handle": + pass + cdef cppclass Simplex_tree_interface_full_featured "Gudhi::Simplex_tree_interface": Simplex_tree() double simplex_filtration(vector[int] simplex) @@ -34,7 +37,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, @@ -43,6 +45,10 @@ 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_filtration(Simplex_tree_simplex_handle f_simplex) + vector[Simplex_tree_simplex_handle].const_iterator get_filtration_iterator_begin() + vector[Simplex_tree_simplex_handle].const_iterator get_filtration_iterator_end() cdef extern from "Persistent_cohomology_interface.h" namespace "Gudhi": cdef cppclass Simplex_tree_persistence_interface "Gudhi::Persistent_cohomology_interface>": diff --git a/src/python/gudhi/simplex_tree.pyx b/src/python/gudhi/simplex_tree.pyx index b18627c4..478139de 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 @@ -214,15 +215,14 @@ cdef class SimplexTree: :returns: The simplices sorted by increasing filtration values. :rtype: list of 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 True: + yield(self.get_ptr().get_simplex_filtration(dereference(it))) + preincrement(it) + if it == end: + raise StopIteration def get_skeleton(self, dimension): """This function returns the (simplices of the) skeleton of a maximum diff --git a/src/python/include/Simplex_tree_interface.h b/src/python/include/Simplex_tree_interface.h index 06f31341..843966cd 100644 --- a/src/python/include/Simplex_tree_interface.h +++ b/src/python/include/Simplex_tree_interface.h @@ -33,7 +33,8 @@ class Simplex_tree_interface : public Simplex_tree { using Simplex_handle = typename Base::Simplex_handle; using Insertion_result = typename std::pair; using Simplex = std::vector; - using Filtered_simplices = std::vector>; + using Filtered_simplex = std::pair; + using Filtered_simplices = std::vector; public: bool find_simplex(const Simplex& vh) { @@ -82,17 +83,12 @@ class Simplex_tree_interface : public Simplex_tree { 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))); + Filtered_simplex get_simplex_filtration(Simplex_handle f_simplex) { + Simplex simplex; + for (auto vertex : Base::simplex_vertex_range(f_simplex)) { + simplex.insert(simplex.begin(), vertex); } - return filtrations; + return std::make_pair(simplex, Base::filtration(f_simplex)); } Filtered_simplices get_skeleton(int dimension) { @@ -135,6 +131,16 @@ class Simplex_tree_interface : public Simplex_tree { Base::initialize_filtration(); pcoh = new Gudhi::Persistent_cohomology_interface(*this); } + + // Iterator over the simplex tree + typename std::vector::const_iterator get_filtration_iterator_begin() { + Base::initialize_filtration(); + return Base::filtration_simplex_range().begin(); + } + + typename std::vector::const_iterator get_filtration_iterator_end() { + return Base::filtration_simplex_range().end(); + } }; } // namespace Gudhi -- cgit v1.2.3 From ef2c5b53e88321f07ad93496f00dde16dc20f018 Mon Sep 17 00:00:00 2001 From: ROUVREAU Vincent Date: Tue, 28 Jan 2020 11:05:39 +0100 Subject: Code review: rename get_simplex_filtration with get_simplex_and_filtration. Remove exception raise. Fix failed tests. Reword documentation --- .../example/alpha_complex_from_points_example.py | 5 +- .../example/rips_complex_from_points_example.py | 5 +- src/python/example/simplex_tree_example.py | 5 +- src/python/gudhi/simplex_tree.pxd | 2 +- src/python/gudhi/simplex_tree.pyx | 10 +-- src/python/include/Simplex_tree_interface.h | 6 +- src/python/test/test_alpha_complex.py | 50 ++++++------ src/python/test/test_euclidean_witness_complex.py | 46 ++++++----- src/python/test/test_rips_complex.py | 53 +++++++------ src/python/test/test_simplex_tree.py | 90 +++++++++++----------- src/python/test/test_tangential_complex.py | 19 +++-- 11 files changed, 161 insertions(+), 130 deletions(-) (limited to 'src/python') diff --git a/src/python/example/alpha_complex_from_points_example.py b/src/python/example/alpha_complex_from_points_example.py index 844d7a82..465632eb 100755 --- a/src/python/example/alpha_complex_from_points_example.py +++ b/src/python/example/alpha_complex_from_points_example.py @@ -47,7 +47,10 @@ else: print("[4] Not found...") 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..7f20c389 100755 --- a/src/python/example/simplex_tree_example.py +++ b/src/python/example/simplex_tree_example.py @@ -39,7 +39,10 @@ else: print("dimension=", st.dimension()) 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/simplex_tree.pxd b/src/python/gudhi/simplex_tree.pxd index caf3c459..1b0dc881 100644 --- a/src/python/gudhi/simplex_tree.pxd +++ b/src/python/gudhi/simplex_tree.pxd @@ -46,7 +46,7 @@ cdef extern from "Simplex_tree_interface.h" namespace "Gudhi": bool prune_above_filtration(double filtration) bool make_filtration_non_decreasing() # Iterators over Simplex tree - pair[vector[int], double] get_simplex_filtration(Simplex_tree_simplex_handle f_simplex) + pair[vector[int], double] get_simplex_and_filtration(Simplex_tree_simplex_handle f_simplex) vector[Simplex_tree_simplex_handle].const_iterator get_filtration_iterator_begin() vector[Simplex_tree_simplex_handle].const_iterator get_filtration_iterator_end() diff --git a/src/python/gudhi/simplex_tree.pyx b/src/python/gudhi/simplex_tree.pyx index 478139de..22978b6e 100644 --- a/src/python/gudhi/simplex_tree.pyx +++ b/src/python/gudhi/simplex_tree.pyx @@ -209,20 +209,18 @@ cdef class SimplexTree: filtration) def get_filtration(self): - """This function returns a list of all simplices with their given + """This function returns a generator with simplices and their given filtration values. :returns: The simplices sorted by increasing filtration values. - :rtype: list of tuples(simplex, filtration) + :rtype: generator with tuples(simplex, filtration) """ 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 True: - yield(self.get_ptr().get_simplex_filtration(dereference(it))) + while it != end: + yield(self.get_ptr().get_simplex_and_filtration(dereference(it))) preincrement(it) - if it == end: - raise StopIteration def get_skeleton(self, dimension): """This function returns the (simplices of the) skeleton of a maximum diff --git a/src/python/include/Simplex_tree_interface.h b/src/python/include/Simplex_tree_interface.h index 843966cd..c0bbc3d9 100644 --- a/src/python/include/Simplex_tree_interface.h +++ b/src/python/include/Simplex_tree_interface.h @@ -33,8 +33,8 @@ class Simplex_tree_interface : public Simplex_tree { using Simplex_handle = typename Base::Simplex_handle; using Insertion_result = typename std::pair; using Simplex = std::vector; - using Filtered_simplex = std::pair; - using Filtered_simplices = std::vector; + using Simplex_and_filtration = std::pair; + using Filtered_simplices = std::vector; public: bool find_simplex(const Simplex& vh) { @@ -83,7 +83,7 @@ class Simplex_tree_interface : public Simplex_tree { Base::initialize_filtration(); } - Filtered_simplex get_simplex_filtration(Simplex_handle 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); diff --git a/src/python/test/test_alpha_complex.py b/src/python/test/test_alpha_complex.py index 3761fe16..ceead919 100755 --- a/src/python/test/test_alpha_complex.py +++ b/src/python/test/test_alpha_complex.py @@ -40,19 +40,21 @@ def test_infinite_alpha(): assert simplex_tree.num_simplices() == 11 assert simplex_tree.num_vertices() == 4 - assert simplex_tree.get_filtration() == [ - ([0], 0.0), - ([1], 0.0), - ([2], 0.0), - ([3], 0.0), - ([0, 1], 0.25), - ([0, 2], 0.25), - ([1, 3], 0.25), - ([2, 3], 0.25), - ([1, 2], 0.5), - ([0, 1, 2], 0.5), - ([1, 2, 3], 0.5), - ] + filtration_generator = simplex_tree.get_filtration() + assert(next(filtration_generator) == ([0], 0.0)) + assert(next(filtration_generator) == ([1], 0.0)) + assert(next(filtration_generator) == ([2], 0.0)) + assert(next(filtration_generator) == ([3], 0.0)) + assert(next(filtration_generator) == ([0, 1], 0.25)) + assert(next(filtration_generator) == ([0, 2], 0.25)) + assert(next(filtration_generator) == ([1, 3], 0.25)) + assert(next(filtration_generator) == ([2, 3], 0.25)) + assert(next(filtration_generator) == ([1, 2], 0.5)) + assert(next(filtration_generator) == ([0, 1, 2], 0.5)) + assert(next(filtration_generator) == ([1, 2, 3], 0.5)) + with pytest.raises(StopIteration): + next(filtration_generator) + assert simplex_tree.get_star([0]) == [ ([0], 0.0), ([0, 1], 0.25), @@ -105,16 +107,18 @@ def test_filtered_alpha(): else: assert False - assert simplex_tree.get_filtration() == [ - ([0], 0.0), - ([1], 0.0), - ([2], 0.0), - ([3], 0.0), - ([0, 1], 0.25), - ([0, 2], 0.25), - ([1, 3], 0.25), - ([2, 3], 0.25), - ] + filtration_generator = simplex_tree.get_filtration() + assert(next(filtration_generator) == ([0], 0.0)) + assert(next(filtration_generator) == ([1], 0.0)) + assert(next(filtration_generator) == ([2], 0.0)) + assert(next(filtration_generator) == ([3], 0.0)) + assert(next(filtration_generator) == ([0, 1], 0.25)) + assert(next(filtration_generator) == ([0, 2], 0.25)) + assert(next(filtration_generator) == ([1, 3], 0.25)) + assert(next(filtration_generator) == ([2, 3], 0.25)) + with pytest.raises(StopIteration): + next(filtration_generator) + assert simplex_tree.get_star([0]) == [([0], 0.0), ([0, 1], 0.25), ([0, 2], 0.25)] assert simplex_tree.get_cofaces([0], 1) == [([0, 1], 0.25), ([0, 2], 0.25)] diff --git a/src/python/test/test_euclidean_witness_complex.py b/src/python/test/test_euclidean_witness_complex.py index c18d2484..16ff1ef4 100755 --- a/src/python/test/test_euclidean_witness_complex.py +++ b/src/python/test/test_euclidean_witness_complex.py @@ -9,6 +9,7 @@ """ import gudhi +import pytest __author__ = "Vincent Rouvreau" __copyright__ = "Copyright (C) 2016 Inria" @@ -40,15 +41,16 @@ 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() == [ - ([0], 0.0), - ([1], 0.0), - ([0, 1], 0.0), - ([2], 0.0), - ([0, 2], 0.0), - ([1, 2], 0.0), - ([0, 1, 2], 0.0), - ] + filtration_generator = simplex_tree.get_filtration() + assert(next(filtration_generator) == ([0], 0.0)) + assert(next(filtration_generator) == ([1], 0.0)) + assert(next(filtration_generator) == ([0, 1], 0.0)) + assert(next(filtration_generator) == ([2], 0.0)) + assert(next(filtration_generator) == ([0, 2], 0.0)) + assert(next(filtration_generator) == ([1, 2], 0.0)) + assert(next(filtration_generator) == ([0, 1, 2], 0.0)) + with pytest.raises(StopIteration): + next(filtration_generator) def test_empty_euclidean_strong_witness_complex(): @@ -78,18 +80,24 @@ 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)] + filtration_generator = simplex_tree.get_filtration() + assert(next(filtration_generator) == ([0], 0.0)) + assert(next(filtration_generator) == ([1], 0.0)) + assert(next(filtration_generator) == ([2], 0.0)) + with pytest.raises(StopIteration): + next(filtration_generator) simplex_tree = euclidean_strong_witness_complex.create_simplex_tree( max_alpha_square=100.0 ) - assert simplex_tree.get_filtration() == [ - ([0], 0.0), - ([1], 0.0), - ([2], 0.0), - ([1, 2], 15.0), - ([0, 2], 34.0), - ([0, 1], 37.0), - ([0, 1, 2], 37.0), - ] + filtration_generator = simplex_tree.get_filtration() + assert(next(filtration_generator) == ([0], 0.0)) + assert(next(filtration_generator) == ([1], 0.0)) + assert(next(filtration_generator) == ([2], 0.0)) + assert(next(filtration_generator) == ([1, 2], 15.0)) + assert(next(filtration_generator) == ([0, 2], 34.0)) + assert(next(filtration_generator) == ([0, 1], 37.0)) + assert(next(filtration_generator) == ([0, 1, 2], 37.0)) + with pytest.raises(StopIteration): + next(filtration_generator) diff --git a/src/python/test/test_rips_complex.py b/src/python/test/test_rips_complex.py index b02a68e1..bd31c47c 100755 --- a/src/python/test/test_rips_complex.py +++ b/src/python/test/test_rips_complex.py @@ -10,6 +10,7 @@ from gudhi import RipsComplex from math import sqrt +import pytest __author__ = "Vincent Rouvreau" __copyright__ = "Copyright (C) 2016 Inria" @@ -32,18 +33,20 @@ def test_rips_from_points(): assert simplex_tree.num_simplices() == 10 assert simplex_tree.num_vertices() == 4 - assert simplex_tree.get_filtration() == [ - ([0], 0.0), - ([1], 0.0), - ([2], 0.0), - ([3], 0.0), - ([0, 1], 1.0), - ([0, 2], 1.0), - ([1, 3], 1.0), - ([2, 3], 1.0), - ([1, 2], 1.4142135623730951), - ([0, 3], 1.4142135623730951), - ] + filtration_generator = simplex_tree.get_filtration() + assert(next(filtration_generator) == ([0], 0.0)) + assert(next(filtration_generator) == ([1], 0.0)) + assert(next(filtration_generator) == ([2], 0.0)) + assert(next(filtration_generator) == ([3], 0.0)) + assert(next(filtration_generator) == ([0, 1], 1.0)) + assert(next(filtration_generator) == ([0, 2], 1.0)) + assert(next(filtration_generator) == ([1, 3], 1.0)) + assert(next(filtration_generator) == ([2, 3], 1.0)) + assert(next(filtration_generator) == ([1, 2], 1.4142135623730951)) + assert(next(filtration_generator) == ([0, 3], 1.4142135623730951)) + with pytest.raises(StopIteration): + next(filtration_generator) + assert simplex_tree.get_star([0]) == [ ([0], 0.0), ([0, 1], 1.0), @@ -95,18 +98,20 @@ def test_rips_from_distance_matrix(): assert simplex_tree.num_simplices() == 10 assert simplex_tree.num_vertices() == 4 - assert simplex_tree.get_filtration() == [ - ([0], 0.0), - ([1], 0.0), - ([2], 0.0), - ([3], 0.0), - ([0, 1], 1.0), - ([0, 2], 1.0), - ([1, 3], 1.0), - ([2, 3], 1.0), - ([1, 2], 1.4142135623730951), - ([0, 3], 1.4142135623730951), - ] + filtration_generator = simplex_tree.get_filtration() + assert(next(filtration_generator) == ([0], 0.0)) + assert(next(filtration_generator) == ([1], 0.0)) + assert(next(filtration_generator) == ([2], 0.0)) + assert(next(filtration_generator) == ([3], 0.0)) + assert(next(filtration_generator) == ([0, 1], 1.0)) + assert(next(filtration_generator) == ([0, 2], 1.0)) + assert(next(filtration_generator) == ([1, 3], 1.0)) + assert(next(filtration_generator) == ([2, 3], 1.0)) + assert(next(filtration_generator) == ([1, 2], 1.4142135623730951)) + assert(next(filtration_generator) == ([0, 3], 1.4142135623730951)) + with pytest.raises(StopIteration): + next(filtration_generator) + 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..0f3db7ac 100755 --- a/src/python/test/test_simplex_tree.py +++ b/src/python/test/test_simplex_tree.py @@ -9,6 +9,7 @@ """ from gudhi import SimplexTree +import pytest __author__ = "Vincent Rouvreau" __copyright__ = "Copyright (C) 2016 Inria" @@ -126,55 +127,58 @@ def test_expansion(): assert st.num_vertices() == 7 assert st.num_simplices() == 17 - assert st.get_filtration() == [ - ([2], 0.1), - ([3], 0.1), - ([2, 3], 0.1), - ([0], 0.2), - ([0, 2], 0.2), - ([1], 0.3), - ([0, 1], 0.3), - ([1, 3], 0.4), - ([1, 2], 0.5), - ([5], 0.6), - ([6], 0.6), - ([5, 6], 0.6), - ([4], 0.7), - ([2, 4], 0.7), - ([0, 3], 0.8), - ([4, 6], 0.9), - ([3, 6], 1.0), - ] + + filtration_generator = st.get_filtration() + assert(next(filtration_generator) == ([2], 0.1)) + assert(next(filtration_generator) == ([3], 0.1)) + assert(next(filtration_generator) == ([2, 3], 0.1)) + assert(next(filtration_generator) == ([0], 0.2)) + assert(next(filtration_generator) == ([0, 2], 0.2)) + assert(next(filtration_generator) == ([1], 0.3)) + assert(next(filtration_generator) == ([0, 1], 0.3)) + assert(next(filtration_generator) == ([1, 3], 0.4)) + assert(next(filtration_generator) == ([1, 2], 0.5)) + assert(next(filtration_generator) == ([5], 0.6)) + assert(next(filtration_generator) == ([6], 0.6)) + assert(next(filtration_generator) == ([5, 6], 0.6)) + assert(next(filtration_generator) == ([4], 0.7)) + assert(next(filtration_generator) == ([2, 4], 0.7)) + assert(next(filtration_generator) == ([0, 3], 0.8)) + assert(next(filtration_generator) == ([4, 6], 0.9)) + assert(next(filtration_generator) == ([3, 6], 1.0)) + with pytest.raises(StopIteration): + next(filtration_generator) st.expansion(3) assert st.num_vertices() == 7 assert st.num_simplices() == 22 st.initialize_filtration() - assert st.get_filtration() == [ - ([2], 0.1), - ([3], 0.1), - ([2, 3], 0.1), - ([0], 0.2), - ([0, 2], 0.2), - ([1], 0.3), - ([0, 1], 0.3), - ([1, 3], 0.4), - ([1, 2], 0.5), - ([0, 1, 2], 0.5), - ([1, 2, 3], 0.5), - ([5], 0.6), - ([6], 0.6), - ([5, 6], 0.6), - ([4], 0.7), - ([2, 4], 0.7), - ([0, 3], 0.8), - ([0, 1, 3], 0.8), - ([0, 2, 3], 0.8), - ([0, 1, 2, 3], 0.8), - ([4, 6], 0.9), - ([3, 6], 1.0), - ] + filtration_generator = st.get_filtration() + assert(next(filtration_generator) == ([2], 0.1)) + assert(next(filtration_generator) == ([3], 0.1)) + assert(next(filtration_generator) == ([2, 3], 0.1)) + assert(next(filtration_generator) == ([0], 0.2)) + assert(next(filtration_generator) == ([0, 2], 0.2)) + assert(next(filtration_generator) == ([1], 0.3)) + assert(next(filtration_generator) == ([0, 1], 0.3)) + assert(next(filtration_generator) == ([1, 3], 0.4)) + assert(next(filtration_generator) == ([1, 2], 0.5)) + assert(next(filtration_generator) == ([0, 1, 2], 0.5)) + assert(next(filtration_generator) == ([1, 2, 3], 0.5)) + assert(next(filtration_generator) == ([5], 0.6)) + assert(next(filtration_generator) == ([6], 0.6)) + assert(next(filtration_generator) == ([5, 6], 0.6)) + assert(next(filtration_generator) == ([4], 0.7)) + assert(next(filtration_generator) == ([2, 4], 0.7)) + assert(next(filtration_generator) == ([0, 3], 0.8)) + assert(next(filtration_generator) == ([0, 1, 3], 0.8)) + assert(next(filtration_generator) == ([0, 2, 3], 0.8)) + assert(next(filtration_generator) == ([0, 1, 2, 3], 0.8)) + assert(next(filtration_generator) == ([4, 6], 0.9)) + assert(next(filtration_generator) == ([3, 6], 1.0)) + with pytest.raises(StopIteration): + next(filtration_generator) def test_automatic_dimension(): diff --git a/src/python/test/test_tangential_complex.py b/src/python/test/test_tangential_complex.py index e650e99c..90e2c75b 100755 --- a/src/python/test/test_tangential_complex.py +++ b/src/python/test/test_tangential_complex.py @@ -9,6 +9,7 @@ """ from gudhi import TangentialComplex, SimplexTree +import pytest __author__ = "Vincent Rouvreau" __copyright__ = "Copyright (C) 2016 Inria" @@ -37,14 +38,16 @@ def test_tangential(): assert st.num_simplices() == 6 assert st.num_vertices() == 4 - assert st.get_filtration() == [ - ([0], 0.0), - ([1], 0.0), - ([2], 0.0), - ([0, 2], 0.0), - ([3], 0.0), - ([1, 3], 0.0), - ] + filtration_generator = st.get_filtration() + assert(next(filtration_generator) == ([0], 0.0)) + assert(next(filtration_generator) == ([1], 0.0)) + assert(next(filtration_generator) == ([2], 0.0)) + assert(next(filtration_generator) == ([0, 2], 0.0)) + assert(next(filtration_generator) == ([3], 0.0)) + assert(next(filtration_generator) == ([1, 3], 0.0)) + with pytest.raises(StopIteration): + next(filtration_generator) + assert st.get_cofaces([0], 1) == [([0, 2], 0.0)] assert point_list[0] == tc.get_point(0) -- cgit v1.2.3 From 68b6e3f3d641cd4a1e86f08bff96e417cc17ac59 Mon Sep 17 00:00:00 2001 From: takeshimeonerespect Date: Fri, 31 Jan 2020 08:08:43 +0100 Subject: timedelay added on fork --- src/python/CMakeLists.txt | 5 +++ src/python/doc/point_cloud.rst | 7 ++++ src/python/gudhi/point_cloud/timedelay.py | 56 +++++++++++++++++++++++++++++++ src/python/test/test_point_cloud.py | 35 +++++++++++++++++++ 4 files changed, 103 insertions(+) create mode 100644 src/python/gudhi/point_cloud/timedelay.py create mode 100755 src/python/test/test_point_cloud.py (limited to 'src/python') diff --git a/src/python/CMakeLists.txt b/src/python/CMakeLists.txt index b558d4c4..b23ec8a9 100644 --- a/src/python/CMakeLists.txt +++ b/src/python/CMakeLists.txt @@ -52,6 +52,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}") @@ -221,6 +222,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 @@ -399,6 +401,9 @@ endif(CGAL_FOUND) add_gudhi_py_test(test_representations) endif() + # Point cloud + add_gudhi_py_test(test_point_cloud) + # Documentation generation is available through sphinx - requires all modules if(SPHINX_PATH) if(MATPLOTLIB_FOUND) diff --git a/src/python/doc/point_cloud.rst b/src/python/doc/point_cloud.rst index d668428a..55c74ff3 100644 --- a/src/python/doc/point_cloud.rst +++ b/src/python/doc/point_cloud.rst @@ -20,3 +20,10 @@ Subsampling :members: :special-members: :show-inheritance: + +TimeDelayEmbedding +------------------ + +.. autoclass:: gudhi.point_cloud.timedelay.TimeDelayEmbedding + :members: + diff --git a/src/python/gudhi/point_cloud/timedelay.py b/src/python/gudhi/point_cloud/timedelay.py new file mode 100644 index 00000000..5c7ba542 --- /dev/null +++ b/src/python/gudhi/point_cloud/timedelay.py @@ -0,0 +1,56 @@ +import numpy as np + +class TimeDelayEmbedding: + """Point cloud transformation class. + + Embeds time-series data in the R^d according to Takens' Embedding 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. + + """ + def __init__(self, dim=3, delay=1, skip=1): + self._dim = dim + self._delay = delay + self._skip = skip + + def __call__(self, *args, **kwargs): + return self.transform(*args, **kwargs) + + def _transform(self, ts): + """Guts of transform method.""" + return ts[ + np.add.outer( + np.arange(0, len(ts)-self._delay*(self._dim-1), self._skip), + np.arange(0, self._dim*self._delay, self._delay)) + ] + + def transform(self, ts): + """Transform method. + + Parameters + ---------- + ts : list[float] or list[list[float]] + A single or multiple time-series data. + + Returns + ------- + point clouds : list[list[float, float, float]] or list[list[list[float, float, float]]] + Makes point cloud every a single time-series data. + """ + ndts = np.array(ts) + if ndts.ndim == 1: + # for single. + return self._transform(ndts).tolist() + else: + # for multiple. + return np.apply_along_axis(self._transform, 1, ndts).tolist() diff --git a/src/python/test/test_point_cloud.py b/src/python/test/test_point_cloud.py new file mode 100755 index 00000000..2ee0c1fb --- /dev/null +++ b/src/python/test/test_point_cloud.py @@ -0,0 +1,35 @@ +from gudhi.point_cloud.timedelay import TimeDelayEmbedding + +def test_normal(): + # Sample array + ts = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10] + # Normal case. + prep = TimeDelayEmbedding() + attractor = prep(ts) + assert (attractor[0] == [1, 2, 3]) + assert (attractor[1] == [2, 3, 4]) + assert (attractor[2] == [3, 4, 5]) + assert (attractor[3] == [4, 5, 6]) + assert (attractor[4] == [5, 6, 7]) + assert (attractor[5] == [6, 7, 8]) + assert (attractor[6] == [7, 8, 9]) + assert (attractor[7] == [8, 9, 10]) + # Delay = 3 + prep = TimeDelayEmbedding(delay=3) + attractor = prep(ts) + assert (attractor[0] == [1, 4, 7]) + assert (attractor[1] == [2, 5, 8]) + assert (attractor[2] == [3, 6, 9]) + assert (attractor[3] == [4, 7, 10]) + # Skip = 3 + prep = TimeDelayEmbedding(skip=3) + attractor = prep(ts) + assert (attractor[0] == [1, 2, 3]) + assert (attractor[1] == [4, 5, 6]) + assert (attractor[2] == [7, 8, 9]) + # Delay = 2 / Skip = 2 + prep = TimeDelayEmbedding(delay=2, skip=2) + attractor = prep(ts) + assert (attractor[0] == [1, 3, 5]) + assert (attractor[1] == [3, 5, 7]) + assert (attractor[2] == [5, 7, 9]) -- cgit v1.2.3 From d6afaa8300daa6204282a7d34df6bea33ea59fd2 Mon Sep 17 00:00:00 2001 From: takeshimeonerespect <58589594+takeshimeonerespect@users.noreply.github.com> Date: Mon, 3 Feb 2020 14:13:52 +0900 Subject: Update timedelay.py --- src/python/gudhi/point_cloud/timedelay.py | 8 ++++++++ 1 file changed, 8 insertions(+) (limited to 'src/python') diff --git a/src/python/gudhi/point_cloud/timedelay.py b/src/python/gudhi/point_cloud/timedelay.py index 5c7ba542..f283916d 100644 --- a/src/python/gudhi/point_cloud/timedelay.py +++ b/src/python/gudhi/point_cloud/timedelay.py @@ -1,3 +1,11 @@ +# 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: -- cgit v1.2.3 From eded147ffffe5b7143cad19ecd134fb7a63991a3 Mon Sep 17 00:00:00 2001 From: takenouchi Date: Tue, 4 Feb 2020 14:08:19 +0900 Subject: change a file name --- src/python/test/test_time_delay.py | 35 +++++++++++++++++++++++++++++++++++ 1 file changed, 35 insertions(+) create mode 100755 src/python/test/test_time_delay.py (limited to 'src/python') diff --git a/src/python/test/test_time_delay.py b/src/python/test/test_time_delay.py new file mode 100755 index 00000000..2ee0c1fb --- /dev/null +++ b/src/python/test/test_time_delay.py @@ -0,0 +1,35 @@ +from gudhi.point_cloud.timedelay import TimeDelayEmbedding + +def test_normal(): + # Sample array + ts = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10] + # Normal case. + prep = TimeDelayEmbedding() + attractor = prep(ts) + assert (attractor[0] == [1, 2, 3]) + assert (attractor[1] == [2, 3, 4]) + assert (attractor[2] == [3, 4, 5]) + assert (attractor[3] == [4, 5, 6]) + assert (attractor[4] == [5, 6, 7]) + assert (attractor[5] == [6, 7, 8]) + assert (attractor[6] == [7, 8, 9]) + assert (attractor[7] == [8, 9, 10]) + # Delay = 3 + prep = TimeDelayEmbedding(delay=3) + attractor = prep(ts) + assert (attractor[0] == [1, 4, 7]) + assert (attractor[1] == [2, 5, 8]) + assert (attractor[2] == [3, 6, 9]) + assert (attractor[3] == [4, 7, 10]) + # Skip = 3 + prep = TimeDelayEmbedding(skip=3) + attractor = prep(ts) + assert (attractor[0] == [1, 2, 3]) + assert (attractor[1] == [4, 5, 6]) + assert (attractor[2] == [7, 8, 9]) + # Delay = 2 / Skip = 2 + prep = TimeDelayEmbedding(delay=2, skip=2) + attractor = prep(ts) + assert (attractor[0] == [1, 3, 5]) + assert (attractor[1] == [3, 5, 7]) + assert (attractor[2] == [5, 7, 9]) -- cgit v1.2.3 From 5ddb724824798fe194a66285e29ea4c5cc2713e2 Mon Sep 17 00:00:00 2001 From: takeshimeonerespect Date: Tue, 4 Feb 2020 14:24:27 +0900 Subject: Delete test_point_cloud.py --- src/python/test/test_point_cloud.py | 35 ----------------------------------- 1 file changed, 35 deletions(-) delete mode 100755 src/python/test/test_point_cloud.py (limited to 'src/python') diff --git a/src/python/test/test_point_cloud.py b/src/python/test/test_point_cloud.py deleted file mode 100755 index 2ee0c1fb..00000000 --- a/src/python/test/test_point_cloud.py +++ /dev/null @@ -1,35 +0,0 @@ -from gudhi.point_cloud.timedelay import TimeDelayEmbedding - -def test_normal(): - # Sample array - ts = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10] - # Normal case. - prep = TimeDelayEmbedding() - attractor = prep(ts) - assert (attractor[0] == [1, 2, 3]) - assert (attractor[1] == [2, 3, 4]) - assert (attractor[2] == [3, 4, 5]) - assert (attractor[3] == [4, 5, 6]) - assert (attractor[4] == [5, 6, 7]) - assert (attractor[5] == [6, 7, 8]) - assert (attractor[6] == [7, 8, 9]) - assert (attractor[7] == [8, 9, 10]) - # Delay = 3 - prep = TimeDelayEmbedding(delay=3) - attractor = prep(ts) - assert (attractor[0] == [1, 4, 7]) - assert (attractor[1] == [2, 5, 8]) - assert (attractor[2] == [3, 6, 9]) - assert (attractor[3] == [4, 7, 10]) - # Skip = 3 - prep = TimeDelayEmbedding(skip=3) - attractor = prep(ts) - assert (attractor[0] == [1, 2, 3]) - assert (attractor[1] == [4, 5, 6]) - assert (attractor[2] == [7, 8, 9]) - # Delay = 2 / Skip = 2 - prep = TimeDelayEmbedding(delay=2, skip=2) - attractor = prep(ts) - assert (attractor[0] == [1, 3, 5]) - assert (attractor[1] == [3, 5, 7]) - assert (attractor[2] == [5, 7, 9]) -- cgit v1.2.3 From 596355344e6205d02110e38a0cb7e0a94e8dbd27 Mon Sep 17 00:00:00 2001 From: takenouchi Date: Thu, 6 Feb 2020 16:00:47 +0900 Subject: modify CMakeLists.txt --- src/python/CMakeLists.txt | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) (limited to 'src/python') diff --git a/src/python/CMakeLists.txt b/src/python/CMakeLists.txt index b23ec8a9..798e2907 100644 --- a/src/python/CMakeLists.txt +++ b/src/python/CMakeLists.txt @@ -401,8 +401,8 @@ endif(CGAL_FOUND) add_gudhi_py_test(test_representations) endif() - # Point cloud - add_gudhi_py_test(test_point_cloud) + # Time Delay + add_gudhi_py_test(test_time_delay) # Documentation generation is available through sphinx - requires all modules if(SPHINX_PATH) -- cgit v1.2.3 From 3253abd27129595f7fcd2be4c2285a93aea98690 Mon Sep 17 00:00:00 2001 From: Vincent Rouvreau <10407034+VincentRouvreau@users.noreply.github.com> Date: Tue, 11 Feb 2020 17:05:08 +0100 Subject: Update src/python/gudhi/simplex_tree.pyx Co-Authored-By: Marc Glisse --- src/python/gudhi/simplex_tree.pyx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) (limited to 'src/python') diff --git a/src/python/gudhi/simplex_tree.pyx b/src/python/gudhi/simplex_tree.pyx index 22978b6e..308b3d2d 100644 --- a/src/python/gudhi/simplex_tree.pyx +++ b/src/python/gudhi/simplex_tree.pyx @@ -219,7 +219,7 @@ cdef class SimplexTree: 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))) + yield self.get_ptr().get_simplex_and_filtration(dereference(it)) preincrement(it) def get_skeleton(self, dimension): -- cgit v1.2.3 From 3ea44646f04648d1a456a0fb9526035101fc17ea Mon Sep 17 00:00:00 2001 From: ROUVREAU Vincent Date: Tue, 11 Feb 2020 17:20:24 +0100 Subject: Code review: non-optimal way to test filtration generator --- src/python/test/test_alpha_complex.py | 49 ++++++------- src/python/test/test_euclidean_witness_complex.py | 45 +++++------- src/python/test/test_rips_complex.py | 50 +++++++------ src/python/test/test_simplex_tree.py | 88 +++++++++++------------ src/python/test/test_tangential_complex.py | 17 +++-- 5 files changed, 117 insertions(+), 132 deletions(-) (limited to 'src/python') diff --git a/src/python/test/test_alpha_complex.py b/src/python/test/test_alpha_complex.py index ceead919..77121302 100755 --- a/src/python/test/test_alpha_complex.py +++ b/src/python/test/test_alpha_complex.py @@ -40,20 +40,19 @@ def test_infinite_alpha(): assert simplex_tree.num_simplices() == 11 assert simplex_tree.num_vertices() == 4 - filtration_generator = simplex_tree.get_filtration() - assert(next(filtration_generator) == ([0], 0.0)) - assert(next(filtration_generator) == ([1], 0.0)) - assert(next(filtration_generator) == ([2], 0.0)) - assert(next(filtration_generator) == ([3], 0.0)) - assert(next(filtration_generator) == ([0, 1], 0.25)) - assert(next(filtration_generator) == ([0, 2], 0.25)) - assert(next(filtration_generator) == ([1, 3], 0.25)) - assert(next(filtration_generator) == ([2, 3], 0.25)) - assert(next(filtration_generator) == ([1, 2], 0.5)) - assert(next(filtration_generator) == ([0, 1, 2], 0.5)) - assert(next(filtration_generator) == ([1, 2, 3], 0.5)) - with pytest.raises(StopIteration): - next(filtration_generator) + assert list(simplex_tree.get_filtration()) == [ + ([0], 0.0), + ([1], 0.0), + ([2], 0.0), + ([3], 0.0), + ([0, 1], 0.25), + ([0, 2], 0.25), + ([1, 3], 0.25), + ([2, 3], 0.25), + ([1, 2], 0.5), + ([0, 1, 2], 0.5), + ([1, 2, 3], 0.5), + ] assert simplex_tree.get_star([0]) == [ ([0], 0.0), @@ -107,18 +106,16 @@ def test_filtered_alpha(): else: assert False - filtration_generator = simplex_tree.get_filtration() - assert(next(filtration_generator) == ([0], 0.0)) - assert(next(filtration_generator) == ([1], 0.0)) - assert(next(filtration_generator) == ([2], 0.0)) - assert(next(filtration_generator) == ([3], 0.0)) - assert(next(filtration_generator) == ([0, 1], 0.25)) - assert(next(filtration_generator) == ([0, 2], 0.25)) - assert(next(filtration_generator) == ([1, 3], 0.25)) - assert(next(filtration_generator) == ([2, 3], 0.25)) - with pytest.raises(StopIteration): - next(filtration_generator) - + assert list(simplex_tree.get_filtration()) == [ + ([0], 0.0), + ([1], 0.0), + ([2], 0.0), + ([3], 0.0), + ([0, 1], 0.25), + ([0, 2], 0.25), + ([1, 3], 0.25), + ([2, 3], 0.25), + ] assert simplex_tree.get_star([0]) == [([0], 0.0), ([0, 1], 0.25), ([0, 2], 0.25)] assert simplex_tree.get_cofaces([0], 1) == [([0, 1], 0.25), ([0, 2], 0.25)] diff --git a/src/python/test/test_euclidean_witness_complex.py b/src/python/test/test_euclidean_witness_complex.py index 16ff1ef4..47196a2a 100755 --- a/src/python/test/test_euclidean_witness_complex.py +++ b/src/python/test/test_euclidean_witness_complex.py @@ -41,16 +41,15 @@ def test_witness_complex(): assert landmarks[1] == euclidean_witness_complex.get_point(1) assert landmarks[2] == euclidean_witness_complex.get_point(2) - filtration_generator = simplex_tree.get_filtration() - assert(next(filtration_generator) == ([0], 0.0)) - assert(next(filtration_generator) == ([1], 0.0)) - assert(next(filtration_generator) == ([0, 1], 0.0)) - assert(next(filtration_generator) == ([2], 0.0)) - assert(next(filtration_generator) == ([0, 2], 0.0)) - assert(next(filtration_generator) == ([1, 2], 0.0)) - assert(next(filtration_generator) == ([0, 1, 2], 0.0)) - with pytest.raises(StopIteration): - next(filtration_generator) + assert list(simplex_tree.get_filtration()) == [ + ([0], 0.0), + ([1], 0.0), + ([0, 1], 0.0), + ([2], 0.0), + ([0, 2], 0.0), + ([1, 2], 0.0), + ([0, 1, 2], 0.0), + ] def test_empty_euclidean_strong_witness_complex(): @@ -80,24 +79,18 @@ 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) - filtration_generator = simplex_tree.get_filtration() - assert(next(filtration_generator) == ([0], 0.0)) - assert(next(filtration_generator) == ([1], 0.0)) - assert(next(filtration_generator) == ([2], 0.0)) - with pytest.raises(StopIteration): - next(filtration_generator) + 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 ) - filtration_generator = simplex_tree.get_filtration() - assert(next(filtration_generator) == ([0], 0.0)) - assert(next(filtration_generator) == ([1], 0.0)) - assert(next(filtration_generator) == ([2], 0.0)) - assert(next(filtration_generator) == ([1, 2], 15.0)) - assert(next(filtration_generator) == ([0, 2], 34.0)) - assert(next(filtration_generator) == ([0, 1], 37.0)) - assert(next(filtration_generator) == ([0, 1, 2], 37.0)) - with pytest.raises(StopIteration): - next(filtration_generator) + assert list(simplex_tree.get_filtration()) == [ + ([0], 0.0), + ([1], 0.0), + ([2], 0.0), + ([1, 2], 15.0), + ([0, 2], 34.0), + ([0, 1], 37.0), + ([0, 1, 2], 37.0), + ] diff --git a/src/python/test/test_rips_complex.py b/src/python/test/test_rips_complex.py index bd31c47c..f5c086cb 100755 --- a/src/python/test/test_rips_complex.py +++ b/src/python/test/test_rips_complex.py @@ -33,19 +33,18 @@ def test_rips_from_points(): assert simplex_tree.num_simplices() == 10 assert simplex_tree.num_vertices() == 4 - filtration_generator = simplex_tree.get_filtration() - assert(next(filtration_generator) == ([0], 0.0)) - assert(next(filtration_generator) == ([1], 0.0)) - assert(next(filtration_generator) == ([2], 0.0)) - assert(next(filtration_generator) == ([3], 0.0)) - assert(next(filtration_generator) == ([0, 1], 1.0)) - assert(next(filtration_generator) == ([0, 2], 1.0)) - assert(next(filtration_generator) == ([1, 3], 1.0)) - assert(next(filtration_generator) == ([2, 3], 1.0)) - assert(next(filtration_generator) == ([1, 2], 1.4142135623730951)) - assert(next(filtration_generator) == ([0, 3], 1.4142135623730951)) - with pytest.raises(StopIteration): - next(filtration_generator) + assert list(simplex_tree.get_filtration()) == [ + ([0], 0.0), + ([1], 0.0), + ([2], 0.0), + ([3], 0.0), + ([0, 1], 1.0), + ([0, 2], 1.0), + ([1, 3], 1.0), + ([2, 3], 1.0), + ([1, 2], 1.4142135623730951), + ([0, 3], 1.4142135623730951), + ] assert simplex_tree.get_star([0]) == [ ([0], 0.0), @@ -98,19 +97,18 @@ def test_rips_from_distance_matrix(): assert simplex_tree.num_simplices() == 10 assert simplex_tree.num_vertices() == 4 - filtration_generator = simplex_tree.get_filtration() - assert(next(filtration_generator) == ([0], 0.0)) - assert(next(filtration_generator) == ([1], 0.0)) - assert(next(filtration_generator) == ([2], 0.0)) - assert(next(filtration_generator) == ([3], 0.0)) - assert(next(filtration_generator) == ([0, 1], 1.0)) - assert(next(filtration_generator) == ([0, 2], 1.0)) - assert(next(filtration_generator) == ([1, 3], 1.0)) - assert(next(filtration_generator) == ([2, 3], 1.0)) - assert(next(filtration_generator) == ([1, 2], 1.4142135623730951)) - assert(next(filtration_generator) == ([0, 3], 1.4142135623730951)) - with pytest.raises(StopIteration): - next(filtration_generator) + assert list(simplex_tree.get_filtration()) == [ + ([0], 0.0), + ([1], 0.0), + ([2], 0.0), + ([3], 0.0), + ([0, 1], 1.0), + ([0, 2], 1.0), + ([1, 3], 1.0), + ([2, 3], 1.0), + ([1, 2], 1.4142135623730951), + ([0, 3], 1.4142135623730951), + ] assert simplex_tree.get_star([0]) == [ ([0], 0.0), diff --git a/src/python/test/test_simplex_tree.py b/src/python/test/test_simplex_tree.py index 0f3db7ac..fa42f2ac 100755 --- a/src/python/test/test_simplex_tree.py +++ b/src/python/test/test_simplex_tree.py @@ -128,57 +128,55 @@ def test_expansion(): assert st.num_vertices() == 7 assert st.num_simplices() == 17 - filtration_generator = st.get_filtration() - assert(next(filtration_generator) == ([2], 0.1)) - assert(next(filtration_generator) == ([3], 0.1)) - assert(next(filtration_generator) == ([2, 3], 0.1)) - assert(next(filtration_generator) == ([0], 0.2)) - assert(next(filtration_generator) == ([0, 2], 0.2)) - assert(next(filtration_generator) == ([1], 0.3)) - assert(next(filtration_generator) == ([0, 1], 0.3)) - assert(next(filtration_generator) == ([1, 3], 0.4)) - assert(next(filtration_generator) == ([1, 2], 0.5)) - assert(next(filtration_generator) == ([5], 0.6)) - assert(next(filtration_generator) == ([6], 0.6)) - assert(next(filtration_generator) == ([5, 6], 0.6)) - assert(next(filtration_generator) == ([4], 0.7)) - assert(next(filtration_generator) == ([2, 4], 0.7)) - assert(next(filtration_generator) == ([0, 3], 0.8)) - assert(next(filtration_generator) == ([4, 6], 0.9)) - assert(next(filtration_generator) == ([3, 6], 1.0)) - with pytest.raises(StopIteration): - next(filtration_generator) + assert list(st.get_filtration()) == [ + ([2], 0.1), + ([3], 0.1), + ([2, 3], 0.1), + ([0], 0.2), + ([0, 2], 0.2), + ([1], 0.3), + ([0, 1], 0.3), + ([1, 3], 0.4), + ([1, 2], 0.5), + ([5], 0.6), + ([6], 0.6), + ([5, 6], 0.6), + ([4], 0.7), + ([2, 4], 0.7), + ([0, 3], 0.8), + ([4, 6], 0.9), + ([3, 6], 1.0), + ] st.expansion(3) assert st.num_vertices() == 7 assert st.num_simplices() == 22 st.initialize_filtration() - filtration_generator = st.get_filtration() - assert(next(filtration_generator) == ([2], 0.1)) - assert(next(filtration_generator) == ([3], 0.1)) - assert(next(filtration_generator) == ([2, 3], 0.1)) - assert(next(filtration_generator) == ([0], 0.2)) - assert(next(filtration_generator) == ([0, 2], 0.2)) - assert(next(filtration_generator) == ([1], 0.3)) - assert(next(filtration_generator) == ([0, 1], 0.3)) - assert(next(filtration_generator) == ([1, 3], 0.4)) - assert(next(filtration_generator) == ([1, 2], 0.5)) - assert(next(filtration_generator) == ([0, 1, 2], 0.5)) - assert(next(filtration_generator) == ([1, 2, 3], 0.5)) - assert(next(filtration_generator) == ([5], 0.6)) - assert(next(filtration_generator) == ([6], 0.6)) - assert(next(filtration_generator) == ([5, 6], 0.6)) - assert(next(filtration_generator) == ([4], 0.7)) - assert(next(filtration_generator) == ([2, 4], 0.7)) - assert(next(filtration_generator) == ([0, 3], 0.8)) - assert(next(filtration_generator) == ([0, 1, 3], 0.8)) - assert(next(filtration_generator) == ([0, 2, 3], 0.8)) - assert(next(filtration_generator) == ([0, 1, 2, 3], 0.8)) - assert(next(filtration_generator) == ([4, 6], 0.9)) - assert(next(filtration_generator) == ([3, 6], 1.0)) - with pytest.raises(StopIteration): - next(filtration_generator) + assert list(st.get_filtration()) == [ + ([2], 0.1), + ([3], 0.1), + ([2, 3], 0.1), + ([0], 0.2), + ([0, 2], 0.2), + ([1], 0.3), + ([0, 1], 0.3), + ([1, 3], 0.4), + ([1, 2], 0.5), + ([0, 1, 2], 0.5), + ([1, 2, 3], 0.5), + ([5], 0.6), + ([6], 0.6), + ([5, 6], 0.6), + ([4], 0.7), + ([2, 4], 0.7), + ([0, 3], 0.8), + ([0, 1, 3], 0.8), + ([0, 2, 3], 0.8), + ([0, 1, 2, 3], 0.8), + ([4, 6], 0.9), + ([3, 6], 1.0), + ] def test_automatic_dimension(): diff --git a/src/python/test/test_tangential_complex.py b/src/python/test/test_tangential_complex.py index 90e2c75b..fc500c45 100755 --- a/src/python/test/test_tangential_complex.py +++ b/src/python/test/test_tangential_complex.py @@ -38,15 +38,14 @@ def test_tangential(): assert st.num_simplices() == 6 assert st.num_vertices() == 4 - filtration_generator = st.get_filtration() - assert(next(filtration_generator) == ([0], 0.0)) - assert(next(filtration_generator) == ([1], 0.0)) - assert(next(filtration_generator) == ([2], 0.0)) - assert(next(filtration_generator) == ([0, 2], 0.0)) - assert(next(filtration_generator) == ([3], 0.0)) - assert(next(filtration_generator) == ([1, 3], 0.0)) - with pytest.raises(StopIteration): - next(filtration_generator) + assert list(st.get_filtration()) == [ + ([0], 0.0), + ([1], 0.0), + ([2], 0.0), + ([0, 2], 0.0), + ([3], 0.0), + ([1, 3], 0.0), + ] assert st.get_cofaces([0], 1) == [([0, 2], 0.0)] -- cgit v1.2.3 From 79de1437cb2fa0ab69465a2f2feabe09a12056eb Mon Sep 17 00:00:00 2001 From: Vincent Rouvreau <10407034+VincentRouvreau@users.noreply.github.com> Date: Tue, 11 Feb 2020 17:51:40 +0100 Subject: Update src/python/include/Simplex_tree_interface.h Co-Authored-By: Marc Glisse --- src/python/include/Simplex_tree_interface.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) (limited to 'src/python') diff --git a/src/python/include/Simplex_tree_interface.h b/src/python/include/Simplex_tree_interface.h index c0bbc3d9..878919cc 100644 --- a/src/python/include/Simplex_tree_interface.h +++ b/src/python/include/Simplex_tree_interface.h @@ -88,7 +88,7 @@ class Simplex_tree_interface : public Simplex_tree { for (auto vertex : Base::simplex_vertex_range(f_simplex)) { simplex.insert(simplex.begin(), vertex); } - return std::make_pair(simplex, Base::filtration(f_simplex)); + return std::make_pair(std::move(simplex), Base::filtration(f_simplex)); } Filtered_simplices get_skeleton(int dimension) { -- cgit v1.2.3 From 1edb818b38ace05b230319227e60838b796ddfc5 Mon Sep 17 00:00:00 2001 From: ROUVREAU Vincent Date: Thu, 13 Feb 2020 11:08:44 +0100 Subject: simplex tree skeleton iterator --- src/python/gudhi/simplex_tree.pxd | 10 +++++++++- src/python/gudhi/simplex_tree.pyx | 15 ++++++--------- src/python/include/Simplex_tree_interface.h | 23 ++++++++++------------- src/python/test/test_simplex_tree.py | 8 ++++---- 4 files changed, 29 insertions(+), 27 deletions(-) (limited to 'src/python') diff --git a/src/python/gudhi/simplex_tree.pxd b/src/python/gudhi/simplex_tree.pxd index 1b0dc881..66c173a6 100644 --- a/src/python/gudhi/simplex_tree.pxd +++ b/src/python/gudhi/simplex_tree.pxd @@ -24,6 +24,13 @@ cdef extern from "Simplex_tree_interface.h" namespace "Gudhi": cdef cppclass Simplex_tree_simplex_handle "Gudhi::Simplex_tree_interface::Simplex_handle": pass + cdef cppclass Simplex_tree_skeleton_iterator "Gudhi::Simplex_tree_interface::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": Simplex_tree() double simplex_filtration(vector[int] simplex) @@ -37,7 +44,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_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) @@ -49,6 +55,8 @@ cdef extern from "Simplex_tree_interface.h" namespace "Gudhi": pair[vector[int], double] get_simplex_and_filtration(Simplex_tree_simplex_handle f_simplex) 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>": diff --git a/src/python/gudhi/simplex_tree.pyx b/src/python/gudhi/simplex_tree.pyx index 308b3d2d..efac2d80 100644 --- a/src/python/gudhi/simplex_tree.pyx +++ b/src/python/gudhi/simplex_tree.pyx @@ -231,15 +231,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(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/include/Simplex_tree_interface.h b/src/python/include/Simplex_tree_interface.h index 878919cc..55d5af97 100644 --- a/src/python/include/Simplex_tree_interface.h +++ b/src/python/include/Simplex_tree_interface.h @@ -35,6 +35,7 @@ class Simplex_tree_interface : public Simplex_tree { using Simplex = std::vector; using Simplex_and_filtration = std::pair; using Filtered_simplices = std::vector; + using Skeleton_simplex_iterator = typename Base::Skeleton_simplex_iterator; public: bool find_simplex(const Simplex& vh) { @@ -91,18 +92,6 @@ class Simplex_tree_interface : public Simplex_tree { return std::make_pair(std::move(simplex), Base::filtration(f_simplex)); } - 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))); - } - return skeletons; - } - Filtered_simplices get_star(const Simplex& simplex) { Filtered_simplices star; for (auto f_simplex : Base::star_simplex_range(Base::find(simplex))) { @@ -134,13 +123,21 @@ class Simplex_tree_interface : public Simplex_tree { // Iterator over the simplex tree typename std::vector::const_iterator get_filtration_iterator_begin() { - Base::initialize_filtration(); + // Base::initialize_filtration(); already performed in filtration_simplex_range return Base::filtration_simplex_range().begin(); } typename std::vector::const_iterator get_filtration_iterator_end() { return Base::filtration_simplex_range().end(); } + + Skeleton_simplex_iterator get_skeleton_iterator_begin(int dimension) { + return Base::skeleton_simplex_range(dimension).begin(); + } + + Skeleton_simplex_iterator get_skeleton_iterator_end(int dimension) { + return Base::skeleton_simplex_range(dimension).end(); + } }; } // namespace Gudhi diff --git a/src/python/test/test_simplex_tree.py b/src/python/test/test_simplex_tree.py index fa42f2ac..eca3807b 100755 --- a/src/python/test/test_simplex_tree.py +++ b/src/python/test/test_simplex_tree.py @@ -56,7 +56,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), @@ -65,7 +65,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), @@ -73,12 +73,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), -- cgit v1.2.3 From a6a4f375822cf3e2ca1866d78472e4350140ddbc Mon Sep 17 00:00:00 2001 From: mtakenouchi Date: Fri, 14 Feb 2020 11:02:56 +0900 Subject: Add __init__.py --- src/python/gudhi/point_cloud/__init__.py | 0 1 file changed, 0 insertions(+), 0 deletions(-) create mode 100644 src/python/gudhi/point_cloud/__init__.py (limited to 'src/python') 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 -- cgit v1.2.3 From 9cc9e1cf3cd9ea42908324d410ef68fa12e8e832 Mon Sep 17 00:00:00 2001 From: mtakenouchi Date: Fri, 14 Feb 2020 11:08:50 +0900 Subject: Update timedelay.py --- src/python/gudhi/point_cloud/timedelay.py | 66 ++++++++++++++++++++++--------- 1 file changed, 48 insertions(+), 18 deletions(-) (limited to 'src/python') diff --git a/src/python/gudhi/point_cloud/timedelay.py b/src/python/gudhi/point_cloud/timedelay.py index f283916d..d899da67 100644 --- a/src/python/gudhi/point_cloud/timedelay.py +++ b/src/python/gudhi/point_cloud/timedelay.py @@ -10,30 +10,55 @@ import numpy as np class TimeDelayEmbedding: """Point cloud transformation class. - Embeds time-series data in the R^d according to Takens' Embedding 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. - + Given delay=3 and skip=2, an point cloud which is obtained by embedding + a single time-series data into R^3 is as follows. + + .. code-block:: none + + time-series = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10] + point clouds = [[1, 4, 7], + [3, 6, 9]] + """ def __init__(self, dim=3, delay=1, skip=1): self._dim = dim self._delay = delay self._skip = skip - def __call__(self, *args, **kwargs): - return self.transform(*args, **kwargs) + def __call__(self, ts): + """Transform method for single time-series data. + Parameters + ---------- + ts : list[float] + A single time-series data. + Returns + ------- + point clouds : list[list[float, float, float]] + Makes point cloud every a single time-series data. + Raises + ------- + TypeError + If the parameter's type does not match the desired type. + """ + ndts = np.array(ts) + if ndts.ndim == 1: + return self._transform(ndts) + else: + raise TypeError("Expects 1-dimensional array.") + def fit(self, ts, y=None): + return self + def _transform(self, ts): """Guts of transform method.""" return ts[ @@ -43,22 +68,27 @@ class TimeDelayEmbedding: ] def transform(self, ts): - """Transform method. - + """Transform method for multiple time-series data. Parameters ---------- - ts : list[float] or list[list[float]] - A single or multiple time-series data. - + ts : list[list[float]] + Multiple time-series data. + Attributes + ---------- + ndts : + The ndts means that all time series need to have exactly + the same size. Returns ------- - point clouds : list[list[float, float, float]] or list[list[list[float, float, float]]] + point clouds : list[list[list[float, float, float]]] Makes point cloud every a single time-series data. + Raises + ------- + TypeError + If the parameter's type does not match the desired type. """ ndts = np.array(ts) - if ndts.ndim == 1: - # for single. - return self._transform(ndts).tolist() + if ndts.ndim == 2: + return np.apply_along_axis(self._transform, 1, ndts) else: - # for multiple. - return np.apply_along_axis(self._transform, 1, ndts).tolist() + raise TypeError("Expects 2-dimensional array.") -- cgit v1.2.3 From 2253fd03bb49aea455309f6d633a6edeb2362d79 Mon Sep 17 00:00:00 2001 From: mtakenouchi Date: Fri, 14 Feb 2020 17:52:07 +0900 Subject: Update test_time_delay.py --- src/python/test/test_time_delay.py | 36 ++++++++++++++++++------------------ 1 file changed, 18 insertions(+), 18 deletions(-) (limited to 'src/python') diff --git a/src/python/test/test_time_delay.py b/src/python/test/test_time_delay.py index 2ee0c1fb..d2ffbf40 100755 --- a/src/python/test/test_time_delay.py +++ b/src/python/test/test_time_delay.py @@ -6,30 +6,30 @@ def test_normal(): # Normal case. prep = TimeDelayEmbedding() attractor = prep(ts) - assert (attractor[0] == [1, 2, 3]) - assert (attractor[1] == [2, 3, 4]) - assert (attractor[2] == [3, 4, 5]) - assert (attractor[3] == [4, 5, 6]) - assert (attractor[4] == [5, 6, 7]) - assert (attractor[5] == [6, 7, 8]) - assert (attractor[6] == [7, 8, 9]) - assert (attractor[7] == [8, 9, 10]) + assert (attractor[0] == np.array[1, 2, 3]) + assert (attractor[1] == np.array[2, 3, 4]) + assert (attractor[2] == np.array[3, 4, 5]) + assert (attractor[3] == np.array[4, 5, 6]) + assert (attractor[4] == np.array[5, 6, 7]) + assert (attractor[5] == np.array[6, 7, 8]) + assert (attractor[6] == np.array[7, 8, 9]) + assert (attractor[7] == np.array[8, 9, 10]) # Delay = 3 prep = TimeDelayEmbedding(delay=3) attractor = prep(ts) - assert (attractor[0] == [1, 4, 7]) - assert (attractor[1] == [2, 5, 8]) - assert (attractor[2] == [3, 6, 9]) - assert (attractor[3] == [4, 7, 10]) + assert (attractor[0] == np.array[1, 4, 7]) + assert (attractor[1] == np.array[2, 5, 8]) + assert (attractor[2] == np.array[3, 6, 9]) + assert (attractor[3] == np.array[4, 7, 10]) # Skip = 3 prep = TimeDelayEmbedding(skip=3) attractor = prep(ts) - assert (attractor[0] == [1, 2, 3]) - assert (attractor[1] == [4, 5, 6]) - assert (attractor[2] == [7, 8, 9]) + assert (attractor[0] == np.array[1, 2, 3]) + assert (attractor[1] == np.array[4, 5, 6]) + assert (attractor[2] == np.array[7, 8, 9]) # Delay = 2 / Skip = 2 prep = TimeDelayEmbedding(delay=2, skip=2) attractor = prep(ts) - assert (attractor[0] == [1, 3, 5]) - assert (attractor[1] == [3, 5, 7]) - assert (attractor[2] == [5, 7, 9]) + assert (attractor[0] == np.array[1, 3, 5]) + assert (attractor[1] == np.array[3, 5, 7]) + assert (attractor[2] == np.array[5, 7, 9]) -- cgit v1.2.3 From f58a4120b70487aede3cb4e81fbb15171e34fa37 Mon Sep 17 00:00:00 2001 From: mtakenouchi Date: Fri, 14 Feb 2020 18:24:18 +0900 Subject: Update test_time_delay.py --- src/python/test/test_time_delay.py | 36 ++++++++++++++++++------------------ 1 file changed, 18 insertions(+), 18 deletions(-) (limited to 'src/python') diff --git a/src/python/test/test_time_delay.py b/src/python/test/test_time_delay.py index d2ffbf40..1cdf56f9 100755 --- a/src/python/test/test_time_delay.py +++ b/src/python/test/test_time_delay.py @@ -6,30 +6,30 @@ def test_normal(): # Normal case. prep = TimeDelayEmbedding() attractor = prep(ts) - assert (attractor[0] == np.array[1, 2, 3]) - assert (attractor[1] == np.array[2, 3, 4]) - assert (attractor[2] == np.array[3, 4, 5]) - assert (attractor[3] == np.array[4, 5, 6]) - assert (attractor[4] == np.array[5, 6, 7]) - assert (attractor[5] == np.array[6, 7, 8]) - assert (attractor[6] == np.array[7, 8, 9]) - assert (attractor[7] == np.array[8, 9, 10]) + assert (attractor[0] == np.array([1, 2, 3])) + assert (attractor[1] == np.array([2, 3, 4])) + assert (attractor[2] == np.array([3, 4, 5])) + assert (attractor[3] == np.array([4, 5, 6])) + assert (attractor[4] == np.array([5, 6, 7])) + assert (attractor[5] == np.array([6, 7, 8])) + assert (attractor[6] == np.array([7, 8, 9])) + assert (attractor[7] == np.array([8, 9, 10])) # Delay = 3 prep = TimeDelayEmbedding(delay=3) attractor = prep(ts) - assert (attractor[0] == np.array[1, 4, 7]) - assert (attractor[1] == np.array[2, 5, 8]) - assert (attractor[2] == np.array[3, 6, 9]) - assert (attractor[3] == np.array[4, 7, 10]) + assert (attractor[0] == np.array([1, 4, 7])) + assert (attractor[1] == np.array([2, 5, 8])) + assert (attractor[2] == np.array([3, 6, 9])) + assert (attractor[3] == np.array([4, 7, 10])) # Skip = 3 prep = TimeDelayEmbedding(skip=3) attractor = prep(ts) - assert (attractor[0] == np.array[1, 2, 3]) - assert (attractor[1] == np.array[4, 5, 6]) - assert (attractor[2] == np.array[7, 8, 9]) + assert (attractor[0] == np.array([1, 2, 3])) + assert (attractor[1] == np.array([4, 5, 6])) + assert (attractor[2] == np.array([7, 8, 9])) # Delay = 2 / Skip = 2 prep = TimeDelayEmbedding(delay=2, skip=2) attractor = prep(ts) - assert (attractor[0] == np.array[1, 3, 5]) - assert (attractor[1] == np.array[3, 5, 7]) - assert (attractor[2] == np.array[5, 7, 9]) + assert (attractor[0] == np.array([1, 3, 5])) + assert (attractor[1] == np.array([3, 5, 7])) + assert (attractor[2] == np.array([5, 7, 9])) -- cgit v1.2.3 From 1c0f48fb26bb2e606dfe0a22e62618357686e2c2 Mon Sep 17 00:00:00 2001 From: mtakenouchi Date: Fri, 14 Feb 2020 18:49:27 +0900 Subject: Update test_time_delay.py --- src/python/test/test_time_delay.py | 1 + 1 file changed, 1 insertion(+) (limited to 'src/python') diff --git a/src/python/test/test_time_delay.py b/src/python/test/test_time_delay.py index 1cdf56f9..3b586ad2 100755 --- a/src/python/test/test_time_delay.py +++ b/src/python/test/test_time_delay.py @@ -1,4 +1,5 @@ from gudhi.point_cloud.timedelay import TimeDelayEmbedding +import numpy as np def test_normal(): # Sample array -- cgit v1.2.3 From 39873c0cf43ca7352dddeab8c1cc6a3fc40a2e58 Mon Sep 17 00:00:00 2001 From: mtakenouchi Date: Fri, 14 Feb 2020 19:08:50 +0900 Subject: Update test_time_delay.py --- src/python/test/test_time_delay.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) (limited to 'src/python') diff --git a/src/python/test/test_time_delay.py b/src/python/test/test_time_delay.py index 3b586ad2..7b6562a5 100755 --- a/src/python/test/test_time_delay.py +++ b/src/python/test/test_time_delay.py @@ -7,7 +7,8 @@ def test_normal(): # Normal case. prep = TimeDelayEmbedding() attractor = prep(ts) - assert (attractor[0] == np.array([1, 2, 3])) + assert (attractor[0] == np.array([1, 2, 3]) + print(attractor[0].all())) assert (attractor[1] == np.array([2, 3, 4])) assert (attractor[2] == np.array([3, 4, 5])) assert (attractor[3] == np.array([4, 5, 6])) -- cgit v1.2.3 From 7c6966ee9821aaeb60d282616445a47071ac1fee Mon Sep 17 00:00:00 2001 From: mtakenouchi Date: Fri, 14 Feb 2020 19:20:25 +0900 Subject: Update test_time_delay.py --- src/python/test/test_time_delay.py | 37 ++++++++++++++++++------------------- 1 file changed, 18 insertions(+), 19 deletions(-) (limited to 'src/python') diff --git a/src/python/test/test_time_delay.py b/src/python/test/test_time_delay.py index 7b6562a5..f652fc88 100755 --- a/src/python/test/test_time_delay.py +++ b/src/python/test/test_time_delay.py @@ -7,31 +7,30 @@ def test_normal(): # Normal case. prep = TimeDelayEmbedding() attractor = prep(ts) - assert (attractor[0] == np.array([1, 2, 3]) - print(attractor[0].all())) - assert (attractor[1] == np.array([2, 3, 4])) - assert (attractor[2] == np.array([3, 4, 5])) - assert (attractor[3] == np.array([4, 5, 6])) - assert (attractor[4] == np.array([5, 6, 7])) - assert (attractor[5] == np.array([6, 7, 8])) - assert (attractor[6] == np.array([7, 8, 9])) - assert (attractor[7] == np.array([8, 9, 10])) + assert (attractor[0].all() == np.array([1, 2, 3])) + assert (attractor[1].all() == np.array([2, 3, 4])) + assert (attractor[2].all() == np.array([3, 4, 5])) + assert (attractor[3].all() == np.array([4, 5, 6])) + assert (attractor[4].all() == np.array([5, 6, 7])) + assert (attractor[5].all() == np.array([6, 7, 8])) + assert (attractor[6].all() == np.array([7, 8, 9])) + assert (attractor[7].all() == np.array([8, 9, 10])) # Delay = 3 prep = TimeDelayEmbedding(delay=3) attractor = prep(ts) - assert (attractor[0] == np.array([1, 4, 7])) - assert (attractor[1] == np.array([2, 5, 8])) - assert (attractor[2] == np.array([3, 6, 9])) - assert (attractor[3] == np.array([4, 7, 10])) + assert (attractor[0].all() == np.array([1, 4, 7])) + assert (attractor[1].all() == np.array([2, 5, 8])) + assert (attractor[2].all() == np.array([3, 6, 9])) + assert (attractor[3].all() == np.array([4, 7, 10])) # Skip = 3 prep = TimeDelayEmbedding(skip=3) attractor = prep(ts) - assert (attractor[0] == np.array([1, 2, 3])) - assert (attractor[1] == np.array([4, 5, 6])) - assert (attractor[2] == np.array([7, 8, 9])) + assert (attractor[0].all() == np.array([1, 2, 3])) + assert (attractor[1].all() == np.array([4, 5, 6])) + assert (attractor[2].all() == np.array([7, 8, 9])) # Delay = 2 / Skip = 2 prep = TimeDelayEmbedding(delay=2, skip=2) attractor = prep(ts) - assert (attractor[0] == np.array([1, 3, 5])) - assert (attractor[1] == np.array([3, 5, 7])) - assert (attractor[2] == np.array([5, 7, 9])) + assert (attractor[0].all() == np.array([1, 3, 5])) + assert (attractor[1].all() == np.array([3, 5, 7])) + assert (attractor[2].all() == np.array([5, 7, 9])) -- cgit v1.2.3 From 5023aa0ff30474a96783152844e7fb0ed52e0c98 Mon Sep 17 00:00:00 2001 From: mtakenouchi Date: Fri, 14 Feb 2020 20:25:14 +0900 Subject: Update test_time_delay.py --- src/python/test/test_time_delay.py | 36 ++++++++++++++++++------------------ 1 file changed, 18 insertions(+), 18 deletions(-) (limited to 'src/python') diff --git a/src/python/test/test_time_delay.py b/src/python/test/test_time_delay.py index f652fc88..5464a185 100755 --- a/src/python/test/test_time_delay.py +++ b/src/python/test/test_time_delay.py @@ -7,30 +7,30 @@ def test_normal(): # Normal case. prep = TimeDelayEmbedding() attractor = prep(ts) - assert (attractor[0].all() == np.array([1, 2, 3])) - assert (attractor[1].all() == np.array([2, 3, 4])) - assert (attractor[2].all() == np.array([3, 4, 5])) - assert (attractor[3].all() == np.array([4, 5, 6])) - assert (attractor[4].all() == np.array([5, 6, 7])) - assert (attractor[5].all() == np.array([6, 7, 8])) - assert (attractor[6].all() == np.array([7, 8, 9])) - assert (attractor[7].all() == np.array([8, 9, 10])) + assert (attractor[0] == np.array([1, 2, 3])).all() + assert (attractor[1] == np.array([2, 3, 4])).all() + assert (attractor[2] == np.array([3, 4, 5])).all() + assert (attractor[3] == np.array([4, 5, 6])).all() + assert (attractor[4] == np.array([5, 6, 7])).all() + assert (attractor[5] == np.array([6, 7, 8])).all() + assert (attractor[6] == np.array([7, 8, 9])).all() + assert (attractor[7] == np.array([8, 9, 10])).all() # Delay = 3 prep = TimeDelayEmbedding(delay=3) attractor = prep(ts) - assert (attractor[0].all() == np.array([1, 4, 7])) - assert (attractor[1].all() == np.array([2, 5, 8])) - assert (attractor[2].all() == np.array([3, 6, 9])) - assert (attractor[3].all() == np.array([4, 7, 10])) + assert (attractor[0] == np.array([1, 4, 7])).all() + assert (attractor[1] == np.array([2, 5, 8])).all() + assert (attractor[2] == np.array([3, 6, 9])).all() + assert (attractor[3] == np.array([4, 7, 10])).all() # Skip = 3 prep = TimeDelayEmbedding(skip=3) attractor = prep(ts) - assert (attractor[0].all() == np.array([1, 2, 3])) - assert (attractor[1].all() == np.array([4, 5, 6])) - assert (attractor[2].all() == np.array([7, 8, 9])) + assert (attractor[0] == np.array([1, 2, 3])).all() + assert (attractor[1] == np.array([4, 5, 6])).all() + assert (attractor[2] == np.array([7, 8, 9])).all() # Delay = 2 / Skip = 2 prep = TimeDelayEmbedding(delay=2, skip=2) attractor = prep(ts) - assert (attractor[0].all() == np.array([1, 3, 5])) - assert (attractor[1].all() == np.array([3, 5, 7])) - assert (attractor[2].all() == np.array([5, 7, 9])) + assert (attractor[0] == np.array([1, 3, 5])).all() + assert (attractor[1] == np.array([3, 5, 7])).all() + assert (attractor[2] == np.array([5, 7, 9])).all() -- cgit v1.2.3 From 80d84e5d8f9a24de745d23f7d721ea3e62217ff4 Mon Sep 17 00:00:00 2001 From: mtakenouchi Date: Wed, 19 Feb 2020 12:32:00 +0900 Subject: Update timedelay.py --- src/python/gudhi/point_cloud/timedelay.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) (limited to 'src/python') diff --git a/src/python/gudhi/point_cloud/timedelay.py b/src/python/gudhi/point_cloud/timedelay.py index d899da67..6ad87cdc 100644 --- a/src/python/gudhi/point_cloud/timedelay.py +++ b/src/python/gudhi/point_cloud/timedelay.py @@ -43,7 +43,7 @@ class TimeDelayEmbedding: A single time-series data. Returns ------- - point clouds : list[list[float, float, float]] + point clouds : list of n x 2 numpy arrays Makes point cloud every a single time-series data. Raises ------- @@ -80,7 +80,7 @@ class TimeDelayEmbedding: the same size. Returns ------- - point clouds : list[list[list[float, float, float]]] + point clouds : list of n x 3 numpy arrays Makes point cloud every a single time-series data. Raises ------- -- cgit v1.2.3 From 88964b4ff10798d6d9c3d0a342c004ee6b8b1496 Mon Sep 17 00:00:00 2001 From: mtakenouchi Date: Tue, 25 Feb 2020 13:21:55 +0900 Subject: Update timedelay.py --- src/python/gudhi/point_cloud/timedelay.py | 89 +++++++++++++++---------------- 1 file changed, 44 insertions(+), 45 deletions(-) (limited to 'src/python') diff --git a/src/python/gudhi/point_cloud/timedelay.py b/src/python/gudhi/point_cloud/timedelay.py index 6ad87cdc..d7a1dab7 100644 --- a/src/python/gudhi/point_cloud/timedelay.py +++ b/src/python/gudhi/point_cloud/timedelay.py @@ -8,10 +8,12 @@ import numpy as np + class TimeDelayEmbedding: """Point cloud transformation class. Embeds time-series data in the R^d according to Takens' Embedding Theorem and obtains the coordinates of each point. + Parameters ---------- dim : int, optional (default=3) @@ -20,16 +22,27 @@ class TimeDelayEmbedding: Time-Delay embedding. skip : int, optional (default=1) How often to skip embedded points. - Given delay=3 and skip=2, an point cloud which is obtained by embedding - a single time-series data into R^3 is as follows. - - .. code-block:: none - - time-series = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10] - point clouds = [[1, 4, 7], - [3, 6, 9]] - + + 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 @@ -39,56 +52,42 @@ class TimeDelayEmbedding: """Transform method for single time-series data. Parameters ---------- - ts : list[float] - A single time-series data. + ts : Iterable[float] or Iterable[Iterable[float]] + A single time-series data, with scalar or vector values. + Returns ------- - point clouds : list of n x 2 numpy arrays - Makes point cloud every a single time-series data. - Raises - ------- - TypeError - If the parameter's type does not match the desired type. + point cloud : n x dim numpy arrays + Makes point cloud from a single time-series data. """ - ndts = np.array(ts) - if ndts.ndim == 1: - return self._transform(ndts) - else: - raise TypeError("Expects 1-dimensional array.") + return self._transform(np.array(ts)) def fit(self, ts, y=None): return self def _transform(self, ts): """Guts of transform method.""" - return ts[ - np.add.outer( - np.arange(0, len(ts)-self._delay*(self._dim-1), self._skip), - np.arange(0, self._dim*self._delay, self._delay)) - ] + 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 : list[list[float]] - Multiple time-series data. - Attributes - ---------- - ndts : - The ndts means that all time series need to have exactly - the same size. + 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 3 numpy arrays - Makes point cloud every a single time-series data. - Raises - ------- - TypeError - If the parameter's type does not match the desired type. + point clouds : list of n x dim numpy arrays + Makes point cloud from each time-series data. """ - ndts = np.array(ts) - if ndts.ndim == 2: - return np.apply_along_axis(self._transform, 1, ndts) - else: - raise TypeError("Expects 2-dimensional array.") + return [self._transform(np.array(s)) for s in ts] -- cgit v1.2.3 From 66c96498b994fea1fcaa6877121023410f4209f9 Mon Sep 17 00:00:00 2001 From: mtakenouchi Date: Tue, 25 Feb 2020 13:24:48 +0900 Subject: Update test_time_delay.py --- src/python/test/test_time_delay.py | 51 ++++++++++++++++++++++---------------- 1 file changed, 29 insertions(+), 22 deletions(-) (limited to 'src/python') diff --git a/src/python/test/test_time_delay.py b/src/python/test/test_time_delay.py index 5464a185..1ead9bca 100755 --- a/src/python/test/test_time_delay.py +++ b/src/python/test/test_time_delay.py @@ -1,36 +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() - attractor = prep(ts) - assert (attractor[0] == np.array([1, 2, 3])).all() - assert (attractor[1] == np.array([2, 3, 4])).all() - assert (attractor[2] == np.array([3, 4, 5])).all() - assert (attractor[3] == np.array([4, 5, 6])).all() - assert (attractor[4] == np.array([5, 6, 7])).all() - assert (attractor[5] == np.array([6, 7, 8])).all() - assert (attractor[6] == np.array([7, 8, 9])).all() - assert (attractor[7] == np.array([8, 9, 10])).all() + 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) - attractor = prep(ts) - assert (attractor[0] == np.array([1, 4, 7])).all() - assert (attractor[1] == np.array([2, 5, 8])).all() - assert (attractor[2] == np.array([3, 6, 9])).all() - assert (attractor[3] == np.array([4, 7, 10])).all() + 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) - attractor = prep(ts) - assert (attractor[0] == np.array([1, 2, 3])).all() - assert (attractor[1] == np.array([4, 5, 6])).all() - assert (attractor[2] == np.array([7, 8, 9])).all() + 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) - attractor = prep(ts) - assert (attractor[0] == np.array([1, 3, 5])).all() - assert (attractor[1] == np.array([3, 5, 7])).all() - assert (attractor[2] == np.array([5, 7, 9])).all() + 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() -- cgit v1.2.3 From a74ec878560bbe5fa340b2650ca9c16471b685af Mon Sep 17 00:00:00 2001 From: mtakenouchi Date: Tue, 25 Feb 2020 13:27:03 +0900 Subject: Update point_cloud.rst --- src/python/doc/point_cloud.rst | 1 + 1 file changed, 1 insertion(+) (limited to 'src/python') diff --git a/src/python/doc/point_cloud.rst b/src/python/doc/point_cloud.rst index 55c74ff3..c0d4b303 100644 --- a/src/python/doc/point_cloud.rst +++ b/src/python/doc/point_cloud.rst @@ -26,4 +26,5 @@ TimeDelayEmbedding .. autoclass:: gudhi.point_cloud.timedelay.TimeDelayEmbedding :members: + :special-members: __call__ -- cgit v1.2.3 From f25d0f86fcd4ac9ab2939b2919d7a66df8b21269 Mon Sep 17 00:00:00 2001 From: mtakenouchi Date: Tue, 25 Feb 2020 16:35:41 +0900 Subject: Update timedelay.py --- src/python/gudhi/point_cloud/timedelay.py | 1 + 1 file changed, 1 insertion(+) (limited to 'src/python') diff --git a/src/python/gudhi/point_cloud/timedelay.py b/src/python/gudhi/point_cloud/timedelay.py index d7a1dab7..576f4386 100644 --- a/src/python/gudhi/point_cloud/timedelay.py +++ b/src/python/gudhi/point_cloud/timedelay.py @@ -50,6 +50,7 @@ class TimeDelayEmbedding: def __call__(self, ts): """Transform method for single time-series data. + Parameters ---------- ts : Iterable[float] or Iterable[Iterable[float]] -- cgit v1.2.3 From 2c1edeb7fd241c8718a22618438b482704703b4a Mon Sep 17 00:00:00 2001 From: mtakenouchi Date: Tue, 25 Feb 2020 17:46:28 +0900 Subject: Update timedelay.py --- src/python/gudhi/point_cloud/timedelay.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) (limited to 'src/python') diff --git a/src/python/gudhi/point_cloud/timedelay.py b/src/python/gudhi/point_cloud/timedelay.py index 576f4386..f01df442 100644 --- a/src/python/gudhi/point_cloud/timedelay.py +++ b/src/python/gudhi/point_cloud/timedelay.py @@ -11,8 +11,9 @@ import numpy as np class TimeDelayEmbedding: """Point cloud transformation class. - Embeds time-series data in the R^d according to Takens' Embedding Theorem - and obtains the coordinates of each point. + 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 ---------- -- cgit v1.2.3 From cbb350d81a8c4acadf31b604aaebde209f462e55 Mon Sep 17 00:00:00 2001 From: ROUVREAU Vincent Date: Wed, 26 Feb 2020 09:32:32 +0100 Subject: Code review: remove import pytest leftovers --- src/python/test/test_euclidean_witness_complex.py | 1 - src/python/test/test_rips_complex.py | 1 - src/python/test/test_simplex_tree.py | 1 - src/python/test/test_tangential_complex.py | 1 - 4 files changed, 4 deletions(-) (limited to 'src/python') diff --git a/src/python/test/test_euclidean_witness_complex.py b/src/python/test/test_euclidean_witness_complex.py index 47196a2a..f3664d39 100755 --- a/src/python/test/test_euclidean_witness_complex.py +++ b/src/python/test/test_euclidean_witness_complex.py @@ -9,7 +9,6 @@ """ import gudhi -import pytest __author__ = "Vincent Rouvreau" __copyright__ = "Copyright (C) 2016 Inria" diff --git a/src/python/test/test_rips_complex.py b/src/python/test/test_rips_complex.py index f5c086cb..b86e7498 100755 --- a/src/python/test/test_rips_complex.py +++ b/src/python/test/test_rips_complex.py @@ -10,7 +10,6 @@ from gudhi import RipsComplex from math import sqrt -import pytest __author__ = "Vincent Rouvreau" __copyright__ = "Copyright (C) 2016 Inria" diff --git a/src/python/test/test_simplex_tree.py b/src/python/test/test_simplex_tree.py index eca3807b..04b26e92 100755 --- a/src/python/test/test_simplex_tree.py +++ b/src/python/test/test_simplex_tree.py @@ -9,7 +9,6 @@ """ from gudhi import SimplexTree -import pytest __author__ = "Vincent Rouvreau" __copyright__ = "Copyright (C) 2016 Inria" diff --git a/src/python/test/test_tangential_complex.py b/src/python/test/test_tangential_complex.py index fc500c45..8668a2e0 100755 --- a/src/python/test/test_tangential_complex.py +++ b/src/python/test/test_tangential_complex.py @@ -9,7 +9,6 @@ """ from gudhi import TangentialComplex, SimplexTree -import pytest __author__ = "Vincent Rouvreau" __copyright__ = "Copyright (C) 2016 Inria" -- cgit v1.2.3 From f85742957276cbd15a2724c86cbc7a8279d62ef9 Mon Sep 17 00:00:00 2001 From: ROUVREAU Vincent Date: Wed, 26 Feb 2020 11:11:32 +0100 Subject: Code review: add some comments about range.begin() and range.end() --- src/python/include/Simplex_tree_interface.h | 4 ++++ 1 file changed, 4 insertions(+) (limited to 'src/python') diff --git a/src/python/include/Simplex_tree_interface.h b/src/python/include/Simplex_tree_interface.h index 55d5af97..66ce5afd 100644 --- a/src/python/include/Simplex_tree_interface.h +++ b/src/python/include/Simplex_tree_interface.h @@ -124,18 +124,22 @@ class Simplex_tree_interface : public Simplex_tree { // Iterator over the simplex tree typename std::vector::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::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(); } }; -- cgit v1.2.3 From 8e4f3d151818b78a29d11cdc6ca171947bfd6dd9 Mon Sep 17 00:00:00 2001 From: tlacombe Date: Tue, 3 Mar 2020 15:33:17 +0100 Subject: update wasserstein distance with pot so that it can return optimal matching now! --- src/python/doc/wasserstein_distance_user.rst | 24 ++++++++++ src/python/gudhi/wasserstein.py | 69 ++++++++++++++++++++++------ src/python/test/test_wasserstein_distance.py | 31 +++++++++---- 3 files changed, 102 insertions(+), 22 deletions(-) (limited to 'src/python') diff --git a/src/python/doc/wasserstein_distance_user.rst b/src/python/doc/wasserstein_distance_user.rst index 94b454e2..d3daa318 100644 --- a/src/python/doc/wasserstein_distance_user.rst +++ b/src/python/doc/wasserstein_distance_user.rst @@ -47,3 +47,27 @@ 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 + + diag1 = np.array([[2.7, 3.7],[9.6, 14.],[34.2, 34.974]]) + diag2 = np.array([[2.8, 4.45], [5, 6], [9.5, 14.1]]) + cost, matching = gudhi.wasserstein.wasserstein_distance(diag1, diag2, matching=True, order=1., internal_p=2.) + + message = "Wasserstein distance value = %.2f, optimal matching: %s" %(cost, matching) + print(message) + +The output is: + +.. testoutput:: + + Wasserstein distance value = 2.15, optimal matching: [(0, 0), (1, 2), (2, -1), (-1, 1)] + diff --git a/src/python/gudhi/wasserstein.py b/src/python/gudhi/wasserstein.py index 13102094..ba0f7343 100644 --- a/src/python/gudhi/wasserstein.py +++ b/src/python/gudhi/wasserstein.py @@ -62,14 +62,39 @@ 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 _clean_match(match, n, m): ''' - :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 match: a list of the form [(i,j) ...] + :param n: int, size of the first dgm + :param m: int, size of the second dgm + :return: a modified version of match where indices greater than n, m are replaced by -1, encoding the diagonal. + and (-1, -1) are removed + ''' + new_match = [] + for i,j in match: + if i >= n: + if j < m: + new_match.append((-1, j)) + elif j >= m: + if i < n: + new_match.append((i,-1)) + else: + new_match.append((i,j)) + return new_match + + +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 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... :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 +102,39 @@ 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., [] else: - return _perstot(Y, order, internal_p) + if not matching: + return _perstot(Y, order, internal_p) + else: + return _perstot(Y, order, internal_p), [(-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), [(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[P < 0.5] = 0 # trick to avoid numerical issue, could it be improved? + match = np.argwhere(P) + # Now we turn to -1 points encoding the diagonal + match = _clean_match(match, n, m) + 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/test/test_wasserstein_distance.py b/src/python/test/test_wasserstein_distance.py index 6a6b217b..02a1d2c9 100755 --- a/src/python/test/test_wasserstein_distance.py +++ b/src/python/test/test_wasserstein_distance.py @@ -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 match == [] + match = wasserstein_distance(emptydiag, emptydiag, matching=True, internal_p=np.inf, order=2.24)[1] + assert match == [] + match = wasserstein_distance(emptydiag, diag2, matching=True, internal_p=np.inf, order=2.)[1] + assert match == [(-1, 0), (-1, 1)] + match = wasserstein_distance(diag2, emptydiag, matching=True, internal_p=np.inf, order=2.24)[1] + assert match == [(0, -1), (1, -1)] + match = wasserstein_distance(diag1, diag2, matching=True, internal_p=2., order=2.)[1] + assert 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) -- cgit v1.2.3 From 2141ef8adfee531f3eaf822cf4076b9b010e6f94 Mon Sep 17 00:00:00 2001 From: tlacombe Date: Tue, 3 Mar 2020 16:22:48 +0100 Subject: correction missing arg in test_wasserstein_distance --- src/python/test/test_wasserstein_distance.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) (limited to 'src/python') diff --git a/src/python/test/test_wasserstein_distance.py b/src/python/test/test_wasserstein_distance.py index 02a1d2c9..d0f0323c 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]]) -- cgit v1.2.3 From 570a9b83eb3f714bc52735dae289a5195874bf41 Mon Sep 17 00:00:00 2001 From: tlacombe Date: Thu, 5 Mar 2020 15:40:45 +0100 Subject: completed as... --- src/python/gudhi/wasserstein.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) (limited to 'src/python') diff --git a/src/python/gudhi/wasserstein.py b/src/python/gudhi/wasserstein.py index ba0f7343..aab0cb3c 100644 --- a/src/python/gudhi/wasserstein.py +++ b/src/python/gudhi/wasserstein.py @@ -30,7 +30,9 @@ 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. + 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). ''' Xdiag = _proj_on_diag(X) @@ -88,7 +90,9 @@ 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 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... + :param matching: if True, computes and returns the optimal matching between X and Y, encoded as + a list of tuple [...(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). -- cgit v1.2.3 From d1d25b4ae8d0f778f0e2b3f98449d7d13e466013 Mon Sep 17 00:00:00 2001 From: ROUVREAU Vincent Date: Tue, 10 Mar 2020 09:04:45 +0100 Subject: Fix example - only fails on OSx --- src/python/example/alpha_complex_from_points_example.py | 3 +++ 1 file changed, 3 insertions(+) (limited to 'src/python') diff --git a/src/python/example/alpha_complex_from_points_example.py b/src/python/example/alpha_complex_from_points_example.py index 465632eb..73faf17c 100755 --- a/src/python/example/alpha_complex_from_points_example.py +++ b/src/python/example/alpha_complex_from_points_example.py @@ -46,6 +46,9 @@ 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=") for simplex_with_filtration in simplex_tree.get_filtration(): -- cgit v1.2.3 From 2eca5c75b1fbd7157e2656b875e730dc5f00f373 Mon Sep 17 00:00:00 2001 From: tlacombe Date: Tue, 10 Mar 2020 15:45:45 +0100 Subject: removed P[P < 0.5] thresholding ; as it shouldn't happen anymore. --- src/python/gudhi/wasserstein.py | 1 - 1 file changed, 1 deletion(-) (limited to 'src/python') diff --git a/src/python/gudhi/wasserstein.py b/src/python/gudhi/wasserstein.py index aab0cb3c..e28c63e6 100644 --- a/src/python/gudhi/wasserstein.py +++ b/src/python/gudhi/wasserstein.py @@ -130,7 +130,6 @@ def wasserstein_distance(X, Y, matching=False, order=2., internal_p=2.): if matching: P = ot.emd(a=a,b=b,M=M, numItermax=2000000) ot_cost = np.sum(np.multiply(P,M)) - P[P < 0.5] = 0 # trick to avoid numerical issue, could it be improved? match = np.argwhere(P) # Now we turn to -1 points encoding the diagonal match = _clean_match(match, n, m) -- cgit v1.2.3 From 967ceab26b09ad74e0cff0d84429a766af267f6b Mon Sep 17 00:00:00 2001 From: tlacombe Date: Tue, 10 Mar 2020 16:47:09 +0100 Subject: removed _clean_match and changed matching format, it is now a (n x 2) numpy array --- src/python/gudhi/wasserstein.py | 31 ++++++------------------------- 1 file changed, 6 insertions(+), 25 deletions(-) (limited to 'src/python') diff --git a/src/python/gudhi/wasserstein.py b/src/python/gudhi/wasserstein.py index e28c63e6..9efa946e 100644 --- a/src/python/gudhi/wasserstein.py +++ b/src/python/gudhi/wasserstein.py @@ -64,34 +64,13 @@ def _perstot(X, order, internal_p): return (np.sum(np.linalg.norm(X - Xdiag, ord=internal_p, axis=1)**order))**(1./order) -def _clean_match(match, n, m): - ''' - :param match: a list of the form [(i,j) ...] - :param n: int, size of the first dgm - :param m: int, size of the second dgm - :return: a modified version of match where indices greater than n, m are replaced by -1, encoding the diagonal. - and (-1, -1) are removed - ''' - new_match = [] - for i,j in match: - if i >= n: - if j < m: - new_match.append((-1, j)) - elif j >= m: - if i < n: - new_match.append((i,-1)) - else: - new_match.append((i,j)) - return new_match - - 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 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 list of tuple [...(i,j)...], meaning the i-th point in X is matched to + 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); @@ -114,12 +93,12 @@ def wasserstein_distance(X, Y, matching=False, order=2., internal_p=2.): if not matching: return _perstot(Y, order, internal_p) else: - return _perstot(Y, order, internal_p), [(-1, j) for j in range(m)] + return _perstot(Y, order, internal_p), np.array([[-1, j] for j in range(m)]) elif Y.size == 0: if not matching: return _perstot(X, order, internal_p) else: - return _perstot(X, order, internal_p), [(i, -1) for i in range(n)] + 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.ones(n+1) # weight vector of the input diagram. Uniform here. @@ -130,9 +109,11 @@ def wasserstein_distance(X, Y, matching=False, order=2., internal_p=2.): 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 = _clean_match(match, n, m) + 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. -- cgit v1.2.3 From 4aea5deab6ce4cbb491f4c9c2b7e9f023efbbe01 Mon Sep 17 00:00:00 2001 From: tlacombe Date: Tue, 10 Mar 2020 17:41:38 +0100 Subject: changed output of matching as a (n x 2) array, adapted tests and doc --- src/python/doc/wasserstein_distance_user.rst | 2 +- src/python/gudhi/wasserstein.py | 2 +- src/python/test/test_wasserstein_distance.py | 10 +++++----- 3 files changed, 7 insertions(+), 7 deletions(-) (limited to 'src/python') diff --git a/src/python/doc/wasserstein_distance_user.rst b/src/python/doc/wasserstein_distance_user.rst index d3daa318..9519caa6 100644 --- a/src/python/doc/wasserstein_distance_user.rst +++ b/src/python/doc/wasserstein_distance_user.rst @@ -69,5 +69,5 @@ The output is: .. testoutput:: - Wasserstein distance value = 2.15, optimal matching: [(0, 0), (1, 2), (2, -1), (-1, 1)] + Wasserstein distance value = 2.15, optimal matching: [[0, 0], [1, 2], [2, -1], [-1, 1]] diff --git a/src/python/gudhi/wasserstein.py b/src/python/gudhi/wasserstein.py index 9efa946e..9e4dc7d5 100644 --- a/src/python/gudhi/wasserstein.py +++ b/src/python/gudhi/wasserstein.py @@ -88,7 +88,7 @@ def wasserstein_distance(X, Y, matching=False, order=2., internal_p=2.): if not matching: return 0. else: - return 0., [] + return 0., np.array([]) else: if not matching: return _perstot(Y, order, internal_p) diff --git a/src/python/test/test_wasserstein_distance.py b/src/python/test/test_wasserstein_distance.py index d0f0323c..ca9a4a61 100755 --- a/src/python/test/test_wasserstein_distance.py +++ b/src/python/test/test_wasserstein_distance.py @@ -61,15 +61,15 @@ def _basic_wasserstein(wasserstein_distance, delta, test_infinity=True, test_mat if test_matching: match = wasserstein_distance(emptydiag, emptydiag, matching=True, internal_p=1., order=2)[1] - assert match == [] + assert np.array_equal(match, np.array([])) match = wasserstein_distance(emptydiag, emptydiag, matching=True, internal_p=np.inf, order=2.24)[1] - assert match == [] + assert np.array_equal(match, np.array([])) match = wasserstein_distance(emptydiag, diag2, matching=True, internal_p=np.inf, order=2.)[1] - assert match == [(-1, 0), (-1, 1)] + assert np.array_equal(match , np.array([[-1, 0], [-1, 1]])) match = wasserstein_distance(diag2, emptydiag, matching=True, internal_p=np.inf, order=2.24)[1] - assert match == [(0, -1), (1, -1)] + assert np.array_equal(match , np.array([[0, -1], [1, -1]])) match = wasserstein_distance(diag1, diag2, matching=True, internal_p=2., order=2.)[1] - assert match == [(0, 0), (1, 1), (2, -1)] + assert np.array_equal(match, np.array_equal([[0, 0], [1, 1], [2, -1]])) -- cgit v1.2.3 From fc4e10863d103ee6bc22863f48548fe246a3ddd6 Mon Sep 17 00:00:00 2001 From: tlacombe Date: Tue, 10 Mar 2020 18:03:21 +0100 Subject: correction of typo in the doc --- src/python/gudhi/wasserstein.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) (limited to 'src/python') diff --git a/src/python/gudhi/wasserstein.py b/src/python/gudhi/wasserstein.py index 9e4dc7d5..12337780 100644 --- a/src/python/gudhi/wasserstein.py +++ b/src/python/gudhi/wasserstein.py @@ -30,10 +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]) + 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 proj onto the diagonal. - note also that C[n+1, m+1] = 0 (it costs nothing to move from the diagonal to 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) -- cgit v1.2.3 From 6c369a6aa566dfcb8cdb501d0c39eafb32219669 Mon Sep 17 00:00:00 2001 From: tlacombe Date: Tue, 10 Mar 2020 18:08:15 +0100 Subject: fix typo in test_wasserstein_distance --- src/python/test/test_wasserstein_distance.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) (limited to 'src/python') diff --git a/src/python/test/test_wasserstein_distance.py b/src/python/test/test_wasserstein_distance.py index ca9a4a61..f92208c0 100755 --- a/src/python/test/test_wasserstein_distance.py +++ b/src/python/test/test_wasserstein_distance.py @@ -69,7 +69,7 @@ def _basic_wasserstein(wasserstein_distance, delta, test_infinity=True, test_mat match = wasserstein_distance(diag2, emptydiag, matching=True, internal_p=np.inf, order=2.24)[1] assert np.array_equal(match , np.array([[0, -1], [1, -1]])) match = wasserstein_distance(diag1, diag2, matching=True, internal_p=2., order=2.)[1] - assert np.array_equal(match, np.array_equal([[0, 0], [1, 1], [2, -1]])) + assert np.array_equal(match, np.array([[0, 0], [1, 1], [2, -1]])) -- cgit v1.2.3 From 753290475ab6e95c2de1baad97ee6f755a0ce19a Mon Sep 17 00:00:00 2001 From: Théo Lacombe Date: Tue, 10 Mar 2020 18:25:10 +0100 Subject: Update src/python/gudhi/wasserstein.py Co-Authored-By: Marc Glisse --- src/python/gudhi/wasserstein.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) (limited to 'src/python') diff --git a/src/python/gudhi/wasserstein.py b/src/python/gudhi/wasserstein.py index 12337780..83a682df 100644 --- a/src/python/gudhi/wasserstein.py +++ b/src/python/gudhi/wasserstein.py @@ -32,7 +32,7 @@ def _build_dist_matrix(X, Y, order=2., internal_p=2.): :returns: (n+1) x (m+1) np.array encoding the cost matrix C. 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 proj onto the diagonal. + 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) -- cgit v1.2.3 From c9d6e27495c8927d736d593afb0450360b46ccc9 Mon Sep 17 00:00:00 2001 From: tlacombe Date: Tue, 10 Mar 2020 18:55:19 +0100 Subject: fix indentation in wasserstein --- src/python/gudhi/wasserstein.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) (limited to 'src/python') diff --git a/src/python/gudhi/wasserstein.py b/src/python/gudhi/wasserstein.py index 83a682df..3dd993f9 100644 --- a/src/python/gudhi/wasserstein.py +++ b/src/python/gudhi/wasserstein.py @@ -70,13 +70,13 @@ def wasserstein_distance(X, Y, matching=False, order=2., internal_p=2.): (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. + 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). + 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. + 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) -- cgit v1.2.3 From a17a09a2c58bba79e897d0ba00aada05da556967 Mon Sep 17 00:00:00 2001 From: tlacombe Date: Wed, 11 Mar 2020 10:41:53 +0100 Subject: clean test_wasserstein from useless np.array --- src/python/test/test_wasserstein_distance.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) (limited to 'src/python') diff --git a/src/python/test/test_wasserstein_distance.py b/src/python/test/test_wasserstein_distance.py index f92208c0..0d70e11a 100755 --- a/src/python/test/test_wasserstein_distance.py +++ b/src/python/test/test_wasserstein_distance.py @@ -61,15 +61,15 @@ def _basic_wasserstein(wasserstein_distance, delta, test_infinity=True, test_mat if test_matching: match = wasserstein_distance(emptydiag, emptydiag, matching=True, internal_p=1., order=2)[1] - assert np.array_equal(match, np.array([])) + 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, np.array([])) + assert np.array_equal(match, []) match = wasserstein_distance(emptydiag, diag2, matching=True, internal_p=np.inf, order=2.)[1] - assert np.array_equal(match , np.array([[-1, 0], [-1, 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 , np.array([[0, -1], [1, -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, np.array([[0, 0], [1, 1], [2, -1]])) + assert np.array_equal(match, [[0, 0], [1, 1], [2, -1]]) -- cgit v1.2.3 From 6ed2a97421a223b4ebe31b91f48d779c2209f470 Mon Sep 17 00:00:00 2001 From: ROUVREAU Vincent Date: Mon, 16 Mar 2020 13:38:18 +0100 Subject: Add get_simplices method - contrary to get_filtration method, sort is not performed --- src/Simplex_tree/example/simple_simplex_tree.cpp | 13 +++++++++++-- src/python/example/simplex_tree_example.py | 4 ++++ src/python/gudhi/simplex_tree.pxd | 8 ++++++++ src/python/gudhi/simplex_tree.pyx | 17 ++++++++++++++++- src/python/include/Simplex_tree_interface.h | 11 +++++++++++ src/python/test/test_simplex_tree.py | 12 ++++++++++++ 6 files changed, 62 insertions(+), 3 deletions(-) (limited to 'src/python') diff --git a/src/Simplex_tree/example/simple_simplex_tree.cpp b/src/Simplex_tree/example/simple_simplex_tree.cpp index 4353939f..47ea7e36 100644 --- a/src/Simplex_tree/example/simple_simplex_tree.cpp +++ b/src/Simplex_tree/example/simple_simplex_tree.cpp @@ -166,10 +166,19 @@ int main(int argc, char* const argv[]) { // ++ GENERAL VARIABLE SET std::cout << "********************************************************************\n"; - // Display the Simplex_tree - Can not be done in the middle of 2 inserts std::cout << "* The complex contains " << simplexTree.num_simplices() << " simplices\n"; std::cout << " - dimension " << simplexTree.dimension() << "\n"; - std::cout << "* Iterator on Simplices in the filtration, with [filtration value]:\n"; + std::cout << "* Iterator on simplices, with [filtration value]:\n"; + for (Simplex_tree::Simplex_handle f_simplex : simplexTree.complex_simplex_range()) { + std::cout << " " + << "[" << simplexTree.filtration(f_simplex) << "] "; + for (auto vertex : simplexTree.simplex_vertex_range(f_simplex)) std::cout << "(" << vertex << ")"; + std::cout << std::endl; + } + + std::cout << "********************************************************************\n"; + // Can not be done in the middle of 2 inserts + std::cout << "* Iterator on simplices sorted by filtration values, with [filtration value]:\n"; for (auto f_simplex : simplexTree.filtration_simplex_range()) { std::cout << " " << "[" << simplexTree.filtration(f_simplex) << "] "; diff --git a/src/python/example/simplex_tree_example.py b/src/python/example/simplex_tree_example.py index 7f20c389..34833899 100755 --- a/src/python/example/simplex_tree_example.py +++ b/src/python/example/simplex_tree_example.py @@ -38,6 +38,10 @@ 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=") for simplex_with_filtration in st.get_filtration(): diff --git a/src/python/gudhi/simplex_tree.pxd b/src/python/gudhi/simplex_tree.pxd index 66c173a6..82f155de 100644 --- a/src/python/gudhi/simplex_tree.pxd +++ b/src/python/gudhi/simplex_tree.pxd @@ -24,6 +24,12 @@ cdef extern from "Simplex_tree_interface.h" namespace "Gudhi": cdef cppclass Simplex_tree_simplex_handle "Gudhi::Simplex_tree_interface::Simplex_handle": pass + cdef cppclass Simplex_tree_simplices_iterator "Gudhi::Simplex_tree_interface::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::Skeleton_simplex_iterator": Simplex_tree_skeleton_iterator() Simplex_tree_simplex_handle& operator*() @@ -53,6 +59,8 @@ cdef extern from "Simplex_tree_interface.h" namespace "Gudhi": 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) diff --git a/src/python/gudhi/simplex_tree.pyx b/src/python/gudhi/simplex_tree.pyx index efac2d80..c01cc905 100644 --- a/src/python/gudhi/simplex_tree.pyx +++ b/src/python/gudhi/simplex_tree.pyx @@ -208,10 +208,25 @@ cdef class SimplexTree: return self.get_ptr().insert_simplex_and_subfaces(csimplex, filtration) - def get_filtration(self): + 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: generator with tuples(simplex, filtration) """ diff --git a/src/python/include/Simplex_tree_interface.h b/src/python/include/Simplex_tree_interface.h index 66ce5afd..4a7062d6 100644 --- a/src/python/include/Simplex_tree_interface.h +++ b/src/python/include/Simplex_tree_interface.h @@ -36,6 +36,7 @@ class Simplex_tree_interface : public Simplex_tree { using Simplex_and_filtration = std::pair; using Filtered_simplices = std::vector; 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) { @@ -122,6 +123,16 @@ class Simplex_tree_interface : public Simplex_tree { } // 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::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 diff --git a/src/python/test/test_simplex_tree.py b/src/python/test/test_simplex_tree.py index 04b26e92..f7848379 100755 --- a/src/python/test/test_simplex_tree.py +++ b/src/python/test/test_simplex_tree.py @@ -249,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] -- cgit v1.2.3 From 5c55e976606b4dd020bd4e21c93ae22143ef5348 Mon Sep 17 00:00:00 2001 From: tlacombe Date: Mon, 16 Mar 2020 18:01:16 +0100 Subject: changed doc of matchings for a more explicit (and hopefully sphinx-valid) version --- src/python/doc/wasserstein_distance_user.rst | 29 ++++++++++++++++++++-------- 1 file changed, 21 insertions(+), 8 deletions(-) (limited to 'src/python') diff --git a/src/python/doc/wasserstein_distance_user.rst b/src/python/doc/wasserstein_distance_user.rst index 9519caa6..4c3b53dd 100644 --- a/src/python/doc/wasserstein_distance_user.rst +++ b/src/python/doc/wasserstein_distance_user.rst @@ -58,16 +58,29 @@ An index of -1 represents the diagonal. 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], [5, 6], [9.5, 14.1]]) - cost, matching = gudhi.wasserstein.wasserstein_distance(diag1, diag2, matching=True, order=1., internal_p=2.) - - message = "Wasserstein distance value = %.2f, optimal matching: %s" %(cost, matching) - print(message) + 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(diag1, diag2, matching=True, order=1., internal_p=2.) + + message_cost = "Wasserstein distance value = %.2f" %cost + print(message_cost) + dgm1_to_diagonal = matchings[np.where(matchings[:,0] == -1)][:,1] + dgm2_to_diagonal = matchings[np.where(matchings[:,1] == -1)][:,0] + 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, optimal matching: [[0, 0], [1, 2], [2, -1], [-1, 1]] - + 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 -- cgit v1.2.3 From 66f0b08a8f8d5006f8d29352c169525cc53a22e6 Mon Sep 17 00:00:00 2001 From: tlacombe Date: Mon, 16 Mar 2020 19:11:30 +0100 Subject: changed typo in doc (diag --> dgm), used integer for order and internal p, simplify th ecode --- src/python/doc/wasserstein_distance_user.rst | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) (limited to 'src/python') diff --git a/src/python/doc/wasserstein_distance_user.rst b/src/python/doc/wasserstein_distance_user.rst index 4c3b53dd..f43b2217 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: @@ -60,12 +60,12 @@ An index of -1 represents the diagonal. 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(diag1, diag2, matching=True, order=1., internal_p=2.) + 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[np.where(matchings[:,0] == -1)][:,1] - dgm2_to_diagonal = matchings[np.where(matchings[:,1] == -1)][:,0] + dgm1_to_diagonal = matching[matching[:,0] == -1, 1] + dgm2_to_diagonal = matching[matching[:,1] == -1, 0] off_diagonal_match = np.delete(matchings, np.where(matchings == -1)[0], axis=0) for i,j in off_diagonal_match: -- cgit v1.2.3 From a253c0c4f54a9a148740ed9c20457df0ea43c842 Mon Sep 17 00:00:00 2001 From: tlacombe Date: Mon, 16 Mar 2020 19:36:07 +0100 Subject: correction typo in usr.rst --- src/python/doc/wasserstein_distance_user.rst | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) (limited to 'src/python') diff --git a/src/python/doc/wasserstein_distance_user.rst b/src/python/doc/wasserstein_distance_user.rst index f43b2217..25e51d68 100644 --- a/src/python/doc/wasserstein_distance_user.rst +++ b/src/python/doc/wasserstein_distance_user.rst @@ -64,8 +64,8 @@ An index of -1 represents the diagonal. message_cost = "Wasserstein distance value = %.2f" %cost print(message_cost) - dgm1_to_diagonal = matching[matching[:,0] == -1, 1] - dgm2_to_diagonal = matching[matching[:,1] == -1, 0] + dgm1_to_diagonal = matchings[matchings[:,0] == -1, 1] + dgm2_to_diagonal = matchings[matchings[:,1] == -1, 0] off_diagonal_match = np.delete(matchings, np.where(matchings == -1)[0], axis=0) for i,j in off_diagonal_match: -- cgit v1.2.3 From 60d11e3f06e08b66e49997f389c4dc01b00b793f Mon Sep 17 00:00:00 2001 From: tlacombe Date: Mon, 16 Mar 2020 21:17:38 +0100 Subject: correction of typo in usr.rst --- src/python/doc/wasserstein_distance_user.rst | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) (limited to 'src/python') diff --git a/src/python/doc/wasserstein_distance_user.rst b/src/python/doc/wasserstein_distance_user.rst index 25e51d68..a9b21fa5 100644 --- a/src/python/doc/wasserstein_distance_user.rst +++ b/src/python/doc/wasserstein_distance_user.rst @@ -64,8 +64,8 @@ An index of -1 represents the diagonal. message_cost = "Wasserstein distance value = %.2f" %cost print(message_cost) - dgm1_to_diagonal = matchings[matchings[:,0] == -1, 1] - dgm2_to_diagonal = matchings[matchings[:,1] == -1, 0] + 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: -- cgit v1.2.3