From 436a7fbe36a9de6a969afd5978c3d496773a8690 Mon Sep 17 00:00:00 2001 From: mathieu Date: Thu, 16 Jan 2020 15:46:41 -0500 Subject: added wrapper functions --- src/python/gudhi/cubical_complex.pyx | 29 ++++++++- src/python/gudhi/simplex_tree.pxd | 1 + src/python/gudhi/simplex_tree.pyx | 28 +++++++- .../include/Persistent_cohomology_interface.h | 76 ++++++++++++++++++++++ 4 files changed, 132 insertions(+), 2 deletions(-) (limited to 'src') diff --git a/src/python/gudhi/cubical_complex.pyx b/src/python/gudhi/cubical_complex.pyx index cbeda014..5562e8a7 100644 --- a/src/python/gudhi/cubical_complex.pyx +++ b/src/python/gudhi/cubical_complex.pyx @@ -31,6 +31,7 @@ cdef extern from "Persistent_cohomology_interface.h" namespace "Gudhi": cdef cppclass Cubical_complex_persistence_interface "Gudhi::Persistent_cohomology_interface>": Cubical_complex_persistence_interface(Bitmap_cubical_complex_base_interface * st, bool persistence_dim_max) vector[pair[int, pair[double, double]]] get_persistence(int homology_coeff_field, double min_persistence) + vector[pair[int, pair[pair[double, int], pair[double, int]]]] get_persistence_cubical_generators(int homology_coeff_field, double min_persistence) vector[int] betti_numbers() vector[int] persistent_betti_numbers(double from_value, double to_value) vector[pair[double,double]] intervals_in_dimension(int dimension) @@ -85,7 +86,7 @@ cdef class CubicalComplex: elif ((dimensions is None) and (top_dimensional_cells is None) and (perseus_file != '')): if os.path.isfile(perseus_file): - self.thisptr = new Bitmap_cubical_complex_base_interface(perseus_file.encode('utf-8')) + self.thisptr = new Bitmap_cubical_complex_base_interface(str.encode(perseus_file)) else: print("file " + perseus_file + " not found.") else: @@ -145,6 +146,32 @@ cdef class CubicalComplex: persistence_result = self.pcohptr.get_persistence(homology_coeff_field, min_persistence) return persistence_result + def persistence_generators(self, homology_coeff_field=11, min_persistence=0, persistence_dim_max = False): + """This function returns the persistence of the simplicial complex. + + :param homology_coeff_field: The homology coefficient field. Must be a + prime number. Default value is 11. + :type homology_coeff_field: int. + :param min_persistence: The minimum persistence value to take into + account (strictly greater than min_persistence). Default value is + 0.0. + Sets min_persistence to -1.0 to see all values. + :type min_persistence: float. + :param persistence_dim_max: If true, the persistent homology for the + maximal dimension in the complex is computed. If false, it is + ignored. Default is false. + :type persistence_dim_max: bool + :returns: The persistence of the simplicial complex, together with the corresponding generators, i.e., the positive and negative top-dimensional cells. + :rtype: list of pairs(dimension, pair(index of positive top-dimensional cell, index of negative top-dimensional cell)) + """ + if self.pcohptr != NULL: + del self.pcohptr + self.pcohptr = new Cubical_complex_persistence_interface(self.thisptr, True) + cdef vector[pair[int, pair[pair[double, int], pair[double, int]]]] persistence_result + if self.pcohptr != NULL: + persistence_result = self.pcohptr.get_persistence_cubical_generators(homology_coeff_field, min_persistence) + return persistence_result + def betti_numbers(self): """This function returns the Betti numbers of the complex. diff --git a/src/python/gudhi/simplex_tree.pxd b/src/python/gudhi/simplex_tree.pxd index 1066d44b..9e52a8aa 100644 --- a/src/python/gudhi/simplex_tree.pxd +++ b/src/python/gudhi/simplex_tree.pxd @@ -48,6 +48,7 @@ cdef extern from "Persistent_cohomology_interface.h" namespace "Gudhi": cdef cppclass Simplex_tree_persistence_interface "Gudhi::Persistent_cohomology_interface>": Simplex_tree_persistence_interface(Simplex_tree_interface_full_featured * st, bool persistence_dim_max) vector[pair[int, pair[double, double]]] get_persistence(int homology_coeff_field, double min_persistence) + vector[pair[int, pair[pair[double, vector[int]], pair[double, vector[int]]]]] get_persistence_generators(int homology_coeff_field, double min_persistence) vector[int] betti_numbers() vector[int] persistent_betti_numbers(double from_value, double to_value) vector[pair[double,double]] intervals_in_dimension(int dimension) diff --git a/src/python/gudhi/simplex_tree.pyx b/src/python/gudhi/simplex_tree.pyx index b18627c4..8cc58f8f 100644 --- a/src/python/gudhi/simplex_tree.pyx +++ b/src/python/gudhi/simplex_tree.pyx @@ -412,6 +412,32 @@ cdef class SimplexTree: persistence_result = self.pcohptr.get_persistence(homology_coeff_field, min_persistence) return persistence_result + def persistence_generators(self, homology_coeff_field=11, min_persistence=0, persistence_dim_max = False): + """This function returns the persistence of the simplicial complex. + + :param homology_coeff_field: The homology coefficient field. Must be a + prime number. Default value is 11. + :type homology_coeff_field: int. + :param min_persistence: The minimum persistence value to take into + account (strictly greater than min_persistence). Default value is + 0.0. + Sets min_persistence to -1.0 to see all values. + :type min_persistence: float. + :param persistence_dim_max: If true, the persistent homology for the + maximal dimension in the complex is computed. If false, it is + ignored. Default is false. + :type persistence_dim_max: bool + :returns: The persistence of the simplicial complex, together with the corresponding generators, i.e., the positive and negative simplices. + :rtype: list of pairs(dimension, pair(positive_simplex, negative_simplex)) + """ + if self.pcohptr != NULL: + del self.pcohptr + self.pcohptr = new Simplex_tree_persistence_interface(self.get_ptr(), persistence_dim_max) + cdef vector[pair[int, pair[pair[double, vector[int]], pair[double, vector[int]]]]] persistence_result + if self.pcohptr != NULL: + persistence_result = self.pcohptr.get_persistence_generators(homology_coeff_field, min_persistence) + return persistence_result + def betti_numbers(self): """This function returns the Betti numbers of the simplicial complex. @@ -508,7 +534,7 @@ cdef class SimplexTree: """ if self.pcohptr != NULL: if persistence_file != '': - self.pcohptr.write_output_diagram(persistence_file.encode('utf-8')) + self.pcohptr.write_output_diagram(str.encode(persistence_file)) else: print("persistence_file must be specified") else: diff --git a/src/python/include/Persistent_cohomology_interface.h b/src/python/include/Persistent_cohomology_interface.h index 8c79e6f3..774eb56a 100644 --- a/src/python/include/Persistent_cohomology_interface.h +++ b/src/python/include/Persistent_cohomology_interface.h @@ -73,6 +73,82 @@ persistent_cohomology::Persistent_cohomology>, std::pair>>>> get_persistence_generators(int homology_coeff_field, + double min_persistence) { + persistent_cohomology::Persistent_cohomology::init_coefficients(homology_coeff_field); + persistent_cohomology::Persistent_cohomology::compute_persistent_cohomology(min_persistence); + + // Custom sort and output persistence + cmp_intervals_by_dim_then_length cmp(stptr_); + auto persistent_pairs = persistent_cohomology::Persistent_cohomology::get_persistent_pairs(); + std::sort(std::begin(persistent_pairs), std::end(persistent_pairs), cmp); + + std::vector>, std::pair>>>> persistence; + for (auto pair : persistent_pairs) { + std::vector splx0, splx1; + for (auto vertex : stptr_->simplex_vertex_range(get<0>(pair))){splx0.push_back(vertex);} + if (isfinite(stptr_->filtration(get<1>(pair)))){ for (auto vertex : stptr_->simplex_vertex_range(get<1>(pair))){splx1.push_back(vertex);}} + persistence.push_back(std::make_pair(stptr_->dimension(get<0>(pair)), std::make_pair(std::make_pair(stptr_->filtration(get<0>(pair)), splx0), std::make_pair(stptr_->filtration(get<1>(pair)), splx1)))); + } + return persistence; + } + + void top_dimensional_cofaces(std::vector & cofaces, int splx){ + if (stptr_->dimension(stptr_->simplex(splx)) == stptr_->dimension()){cofaces.push_back(stptr_->simplex(splx));} + else{ for (auto v : stptr_->coboundary_simplex_range(stptr_->simplex(splx))){top_dimensional_cofaces(cofaces, stptr_->key(v));} } + } + + std::vector, std::pair>>> get_persistence_cubical_generators(int homology_coeff_field, + double min_persistence) { + + // Gather all top-dimensional cells and store their simplex handles + std::vector max_splx; for (auto splx : stptr_->filtration_simplex_range()){ if (stptr_->dimension(splx) == stptr_->dimension()) max_splx.push_back(splx); } + // Sort these simplex handles and compute the ordering function + // This function allows to go directly from the simplex handle to the position of the corresponding top-dimensional cell in the input data + std::map order; std::sort(max_splx.begin(), max_splx.end()); for (int i = 0; i < max_splx.size(); i++) order.insert(std::make_pair(max_splx[i], i)); + + persistent_cohomology::Persistent_cohomology::init_coefficients(homology_coeff_field); + persistent_cohomology::Persistent_cohomology::compute_persistent_cohomology(min_persistence); + + // Custom sort and output persistence + cmp_intervals_by_dim_then_length cmp(stptr_); + auto persistent_pairs = persistent_cohomology::Persistent_cohomology::get_persistent_pairs(); + std::sort(std::begin(persistent_pairs), std::end(persistent_pairs), cmp); + + std::vector, std::pair>>> persistence; + for (auto pair : persistent_pairs) { + + int splx0, splx1; + + double f0 = stptr_->filtration(get<0>(pair)); + // Recursively get the top-dimensional cells / cofaces associated to the persistence generator + std::vector faces0; top_dimensional_cofaces(faces0, stptr_->key(get<0>(pair))); + // Find the top-dimensional cell / coface with the same filtration value + int cf; for (int i = 0; i < faces0.size(); i++){ if (stptr_->filtration(faces0[i]) == f0){cf = i; break;}} + // Retrieve the index of the corresponding top-dimensional cell in the input data + splx0 = order[faces0[cf]]; + + if (isfinite(stptr_->filtration(get<1>(pair)))){ + double f1 = stptr_->filtration(get<1>(pair)); + // Recursively get the top-dimensional cells / cofaces associated to the persistence generator + std::vector faces1; top_dimensional_cofaces(faces1, stptr_->key(get<1>(pair))); + // Find the top-dimensional cell / coface with the same filtration value + int cf; for (int i = 0; i < faces0.size(); i++){ if (stptr_->filtration(faces0[i]) == f0){cf = i; break;}} + // Retrieve the index of the corresponding top-dimensional cell in the input data + splx1 = order[faces1[cf]]; + } + + persistence.push_back(std::make_pair(stptr_->dimension(get<0>(pair)), std::make_pair(std::make_pair(stptr_->filtration(get<0>(pair)), splx0), std::make_pair(stptr_->filtration(get<1>(pair)), splx1)))); + } + return persistence; + } + std::vector, std::vector>> persistence_pairs() { auto pairs = persistent_cohomology::Persistent_cohomology::get_persistent_pairs(); -- cgit v1.2.3 From 5694670b3e20f0cb935a751614ef12b6009a60c0 Mon Sep 17 00:00:00 2001 From: mathieu Date: Thu, 16 Jan 2020 15:58:15 -0500 Subject: fix to detect infinite persistence --- src/python/include/Persistent_cohomology_interface.h | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) (limited to 'src') diff --git a/src/python/include/Persistent_cohomology_interface.h b/src/python/include/Persistent_cohomology_interface.h index 774eb56a..acc32b21 100644 --- a/src/python/include/Persistent_cohomology_interface.h +++ b/src/python/include/Persistent_cohomology_interface.h @@ -124,16 +124,15 @@ persistent_cohomology::Persistent_cohomology, std::pair>>> persistence; for (auto pair : persistent_pairs) { - int splx0, splx1; - double f0 = stptr_->filtration(get<0>(pair)); // Recursively get the top-dimensional cells / cofaces associated to the persistence generator std::vector faces0; top_dimensional_cofaces(faces0, stptr_->key(get<0>(pair))); // Find the top-dimensional cell / coface with the same filtration value int cf; for (int i = 0; i < faces0.size(); i++){ if (stptr_->filtration(faces0[i]) == f0){cf = i; break;}} // Retrieve the index of the corresponding top-dimensional cell in the input data - splx0 = order[faces0[cf]]; + int splx0 = order[faces0[cf]]; + int splx1 = -1; if (isfinite(stptr_->filtration(get<1>(pair)))){ double f1 = stptr_->filtration(get<1>(pair)); // Recursively get the top-dimensional cells / cofaces associated to the persistence generator -- cgit v1.2.3 From c89df405c77bb7270db1a7d8f0e49bc22c1b010d Mon Sep 17 00:00:00 2001 From: mathieu Date: Thu, 16 Jan 2020 16:17:38 -0500 Subject: fix typo + coboundary error --- src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex.h | 1 + src/python/include/Persistent_cohomology_interface.h | 6 +++--- 2 files changed, 4 insertions(+), 3 deletions(-) (limited to 'src') diff --git a/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex.h b/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex.h index 37514dee..bf09532e 100644 --- a/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex.h +++ b/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex.h @@ -340,6 +340,7 @@ class Bitmap_cubical_complex : public T { * that provides ranges for the Boundary_simplex_iterator. **/ Boundary_simplex_range boundary_simplex_range(Simplex_handle sh) { return this->get_boundary_of_a_cell(sh); } + Boundary_simplex_range coboundary_simplex_range(Simplex_handle sh) { return this->get_coboundary_of_a_cell(sh); } /** * filtration_simplex_range creates an object of a Filtration_simplex_range class diff --git a/src/python/include/Persistent_cohomology_interface.h b/src/python/include/Persistent_cohomology_interface.h index acc32b21..0ad14477 100644 --- a/src/python/include/Persistent_cohomology_interface.h +++ b/src/python/include/Persistent_cohomology_interface.h @@ -108,7 +108,7 @@ persistent_cohomology::Persistent_cohomology max_splx; for (auto splx : stptr_->filtration_simplex_range()){ if (stptr_->dimension(splx) == stptr_->dimension()) max_splx.push_back(splx); } // Sort these simplex handles and compute the ordering function // This function allows to go directly from the simplex handle to the position of the corresponding top-dimensional cell in the input data - std::map order; std::sort(max_splx.begin(), max_splx.end()); for (int i = 0; i < max_splx.size(); i++) order.insert(std::make_pair(max_splx[i], i)); + std::map order; std::sort(max_splx.begin(), max_splx.end()); for (unsigned int i = 0; i < max_splx.size(); i++) order.insert(std::make_pair(max_splx[i], i)); persistent_cohomology::Persistent_cohomology::init_coefficients(homology_coeff_field); @@ -128,7 +128,7 @@ persistent_cohomology::Persistent_cohomology faces0; top_dimensional_cofaces(faces0, stptr_->key(get<0>(pair))); // Find the top-dimensional cell / coface with the same filtration value - int cf; for (int i = 0; i < faces0.size(); i++){ if (stptr_->filtration(faces0[i]) == f0){cf = i; break;}} + int cf; for (unsigned int i = 0; i < faces0.size(); i++){if (stptr_->filtration(faces0[i]) == f0){cf = i; break;}} // Retrieve the index of the corresponding top-dimensional cell in the input data int splx0 = order[faces0[cf]]; @@ -138,7 +138,7 @@ persistent_cohomology::Persistent_cohomology faces1; top_dimensional_cofaces(faces1, stptr_->key(get<1>(pair))); // Find the top-dimensional cell / coface with the same filtration value - int cf; for (int i = 0; i < faces0.size(); i++){ if (stptr_->filtration(faces0[i]) == f0){cf = i; break;}} + int cf; for (unsigned int i = 0; i < faces1.size(); i++){if (stptr_->filtration(faces1[i]) == f1){cf = i; break;}} // Retrieve the index of the corresponding top-dimensional cell in the input data splx1 = order[faces1[cf]]; } -- cgit v1.2.3 From 19562b27182dcfa6ed262002c2bc8934382f5a53 Mon Sep 17 00:00:00 2001 From: MathieuCarriere Date: Thu, 16 Jan 2020 21:26:02 -0500 Subject: get rid of persistence_generators and modified name for cubical complex --- src/python/gudhi/cubical_complex.pyx | 8 +++--- src/python/gudhi/simplex_tree.pxd | 1 - src/python/gudhi/simplex_tree.pyx | 26 ------------------- .../include/Persistent_cohomology_interface.h | 29 +++------------------- 4 files changed, 8 insertions(+), 56 deletions(-) (limited to 'src') diff --git a/src/python/gudhi/cubical_complex.pyx b/src/python/gudhi/cubical_complex.pyx index 5562e8a7..8ea31486 100644 --- a/src/python/gudhi/cubical_complex.pyx +++ b/src/python/gudhi/cubical_complex.pyx @@ -31,7 +31,7 @@ cdef extern from "Persistent_cohomology_interface.h" namespace "Gudhi": cdef cppclass Cubical_complex_persistence_interface "Gudhi::Persistent_cohomology_interface>": Cubical_complex_persistence_interface(Bitmap_cubical_complex_base_interface * st, bool persistence_dim_max) vector[pair[int, pair[double, double]]] get_persistence(int homology_coeff_field, double min_persistence) - vector[pair[int, pair[pair[double, int], pair[double, int]]]] get_persistence_cubical_generators(int homology_coeff_field, double min_persistence) + vector[pair[int, pair[pair[double, int], pair[double, int]]]] get_cofaces_of_cubical_persistence_pairs(int homology_coeff_field, double min_persistence) vector[int] betti_numbers() vector[int] persistent_betti_numbers(double from_value, double to_value) vector[pair[double,double]] intervals_in_dimension(int dimension) @@ -146,7 +146,7 @@ cdef class CubicalComplex: persistence_result = self.pcohptr.get_persistence(homology_coeff_field, min_persistence) return persistence_result - def persistence_generators(self, homology_coeff_field=11, min_persistence=0, persistence_dim_max = False): + def cofaces_of_cubical_persistence_pairs(self, homology_coeff_field=11, min_persistence=0, persistence_dim_max = False): """This function returns the persistence of the simplicial complex. :param homology_coeff_field: The homology coefficient field. Must be a @@ -161,7 +161,7 @@ cdef class CubicalComplex: maximal dimension in the complex is computed. If false, it is ignored. Default is false. :type persistence_dim_max: bool - :returns: The persistence of the simplicial complex, together with the corresponding generators, i.e., the positive and negative top-dimensional cells. + :returns: The persistence of the simplicial complex, together with the cofaces of the corresponding generators, i.e., the top-dimensional cells/cofaces of the positive and negative simplices. :rtype: list of pairs(dimension, pair(index of positive top-dimensional cell, index of negative top-dimensional cell)) """ if self.pcohptr != NULL: @@ -169,7 +169,7 @@ cdef class CubicalComplex: self.pcohptr = new Cubical_complex_persistence_interface(self.thisptr, True) cdef vector[pair[int, pair[pair[double, int], pair[double, int]]]] persistence_result if self.pcohptr != NULL: - persistence_result = self.pcohptr.get_persistence_cubical_generators(homology_coeff_field, min_persistence) + persistence_result = self.pcohptr.get_cofaces_of_cubical_persistence_pairs(homology_coeff_field, min_persistence) return persistence_result def betti_numbers(self): diff --git a/src/python/gudhi/simplex_tree.pxd b/src/python/gudhi/simplex_tree.pxd index 9e52a8aa..1066d44b 100644 --- a/src/python/gudhi/simplex_tree.pxd +++ b/src/python/gudhi/simplex_tree.pxd @@ -48,7 +48,6 @@ cdef extern from "Persistent_cohomology_interface.h" namespace "Gudhi": cdef cppclass Simplex_tree_persistence_interface "Gudhi::Persistent_cohomology_interface>": Simplex_tree_persistence_interface(Simplex_tree_interface_full_featured * st, bool persistence_dim_max) vector[pair[int, pair[double, double]]] get_persistence(int homology_coeff_field, double min_persistence) - vector[pair[int, pair[pair[double, vector[int]], pair[double, vector[int]]]]] get_persistence_generators(int homology_coeff_field, double min_persistence) vector[int] betti_numbers() vector[int] persistent_betti_numbers(double from_value, double to_value) vector[pair[double,double]] intervals_in_dimension(int dimension) diff --git a/src/python/gudhi/simplex_tree.pyx b/src/python/gudhi/simplex_tree.pyx index 8cc58f8f..85d25492 100644 --- a/src/python/gudhi/simplex_tree.pyx +++ b/src/python/gudhi/simplex_tree.pyx @@ -412,32 +412,6 @@ cdef class SimplexTree: persistence_result = self.pcohptr.get_persistence(homology_coeff_field, min_persistence) return persistence_result - def persistence_generators(self, homology_coeff_field=11, min_persistence=0, persistence_dim_max = False): - """This function returns the persistence of the simplicial complex. - - :param homology_coeff_field: The homology coefficient field. Must be a - prime number. Default value is 11. - :type homology_coeff_field: int. - :param min_persistence: The minimum persistence value to take into - account (strictly greater than min_persistence). Default value is - 0.0. - Sets min_persistence to -1.0 to see all values. - :type min_persistence: float. - :param persistence_dim_max: If true, the persistent homology for the - maximal dimension in the complex is computed. If false, it is - ignored. Default is false. - :type persistence_dim_max: bool - :returns: The persistence of the simplicial complex, together with the corresponding generators, i.e., the positive and negative simplices. - :rtype: list of pairs(dimension, pair(positive_simplex, negative_simplex)) - """ - if self.pcohptr != NULL: - del self.pcohptr - self.pcohptr = new Simplex_tree_persistence_interface(self.get_ptr(), persistence_dim_max) - cdef vector[pair[int, pair[pair[double, vector[int]], pair[double, vector[int]]]]] persistence_result - if self.pcohptr != NULL: - persistence_result = self.pcohptr.get_persistence_generators(homology_coeff_field, min_persistence) - return persistence_result - def betti_numbers(self): """This function returns the Betti numbers of the simplicial complex. diff --git a/src/python/include/Persistent_cohomology_interface.h b/src/python/include/Persistent_cohomology_interface.h index 0ad14477..1a1e716e 100644 --- a/src/python/include/Persistent_cohomology_interface.h +++ b/src/python/include/Persistent_cohomology_interface.h @@ -73,36 +73,15 @@ persistent_cohomology::Persistent_cohomology>, std::pair>>>> get_persistence_generators(int homology_coeff_field, - double min_persistence) { - persistent_cohomology::Persistent_cohomology::init_coefficients(homology_coeff_field); - persistent_cohomology::Persistent_cohomology::compute_persistent_cohomology(min_persistence); - - // Custom sort and output persistence - cmp_intervals_by_dim_then_length cmp(stptr_); - auto persistent_pairs = persistent_cohomology::Persistent_cohomology::get_persistent_pairs(); - std::sort(std::begin(persistent_pairs), std::end(persistent_pairs), cmp); - - std::vector>, std::pair>>>> persistence; - for (auto pair : persistent_pairs) { - std::vector splx0, splx1; - for (auto vertex : stptr_->simplex_vertex_range(get<0>(pair))){splx0.push_back(vertex);} - if (isfinite(stptr_->filtration(get<1>(pair)))){ for (auto vertex : stptr_->simplex_vertex_range(get<1>(pair))){splx1.push_back(vertex);}} - persistence.push_back(std::make_pair(stptr_->dimension(get<0>(pair)), std::make_pair(std::make_pair(stptr_->filtration(get<0>(pair)), splx0), std::make_pair(stptr_->filtration(get<1>(pair)), splx1)))); - } - return persistence; - } - void top_dimensional_cofaces(std::vector & cofaces, int splx){ if (stptr_->dimension(stptr_->simplex(splx)) == stptr_->dimension()){cofaces.push_back(stptr_->simplex(splx));} else{ for (auto v : stptr_->coboundary_simplex_range(stptr_->simplex(splx))){top_dimensional_cofaces(cofaces, stptr_->key(v));} } } - std::vector, std::pair>>> get_persistence_cubical_generators(int homology_coeff_field, - double min_persistence) { + std::vector, std::pair>>> get_cofaces_of_cubical_persistence_pairs(int homology_coeff_field, + double min_persistence) { + + // Warning: this function is meant to be used with CubicalComplex only!! // Gather all top-dimensional cells and store their simplex handles std::vector max_splx; for (auto splx : stptr_->filtration_simplex_range()){ if (stptr_->dimension(splx) == stptr_->dimension()) max_splx.push_back(splx); } -- cgit v1.2.3 From 62e92e64bd97ec0bd26c31e071228f7d7c78b0e5 Mon Sep 17 00:00:00 2001 From: MathieuCarriere Date: Thu, 16 Jan 2020 21:29:55 -0500 Subject: fixed typo for CubicalComplex --- src/python/gudhi/cubical_complex.pyx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) (limited to 'src') diff --git a/src/python/gudhi/cubical_complex.pyx b/src/python/gudhi/cubical_complex.pyx index 8ea31486..bd432834 100644 --- a/src/python/gudhi/cubical_complex.pyx +++ b/src/python/gudhi/cubical_complex.pyx @@ -86,7 +86,7 @@ cdef class CubicalComplex: elif ((dimensions is None) and (top_dimensional_cells is None) and (perseus_file != '')): if os.path.isfile(perseus_file): - self.thisptr = new Bitmap_cubical_complex_base_interface(str.encode(perseus_file)) + self.thisptr = new Bitmap_cubical_complex_base_interface(perseus_file.encode('utf-8')) else: print("file " + perseus_file + " not found.") else: -- cgit v1.2.3 From fe754ca20cf942e2af186f14e5a3d24e23b6c80e Mon Sep 17 00:00:00 2001 From: mathieu Date: Thu, 13 Feb 2020 19:27:40 -0500 Subject: fix Marc's comments --- src/python/gudhi/cubical_complex.pyx | 49 ++++++++-------- .../include/Persistent_cohomology_interface.h | 67 ++++++++++------------ 2 files changed, 55 insertions(+), 61 deletions(-) (limited to 'src') diff --git a/src/python/gudhi/cubical_complex.pyx b/src/python/gudhi/cubical_complex.pyx index bd432834..8cf43539 100644 --- a/src/python/gudhi/cubical_complex.pyx +++ b/src/python/gudhi/cubical_complex.pyx @@ -31,7 +31,7 @@ cdef extern from "Persistent_cohomology_interface.h" namespace "Gudhi": cdef cppclass Cubical_complex_persistence_interface "Gudhi::Persistent_cohomology_interface>": Cubical_complex_persistence_interface(Bitmap_cubical_complex_base_interface * st, bool persistence_dim_max) vector[pair[int, pair[double, double]]] get_persistence(int homology_coeff_field, double min_persistence) - vector[pair[int, pair[pair[double, int], pair[double, int]]]] get_cofaces_of_cubical_persistence_pairs(int homology_coeff_field, double min_persistence) + vector[vector[int]] cofaces_of_cubical_persistence_pairs() vector[int] betti_numbers() vector[int] persistent_betti_numbers(double from_value, double to_value) vector[pair[double,double]] intervals_in_dimension(int dimension) @@ -146,31 +146,32 @@ cdef class CubicalComplex: persistence_result = self.pcohptr.get_persistence(homology_coeff_field, min_persistence) return persistence_result - def cofaces_of_cubical_persistence_pairs(self, homology_coeff_field=11, min_persistence=0, persistence_dim_max = False): - """This function returns the persistence of the simplicial complex. - - :param homology_coeff_field: The homology coefficient field. Must be a - prime number. Default value is 11. - :type homology_coeff_field: int. - :param min_persistence: The minimum persistence value to take into - account (strictly greater than min_persistence). Default value is - 0.0. - Sets min_persistence to -1.0 to see all values. - :type min_persistence: float. - :param persistence_dim_max: If true, the persistent homology for the - maximal dimension in the complex is computed. If false, it is - ignored. Default is false. - :type persistence_dim_max: bool - :returns: The persistence of the simplicial complex, together with the cofaces of the corresponding generators, i.e., the top-dimensional cells/cofaces of the positive and negative simplices. - :rtype: list of pairs(dimension, pair(index of positive top-dimensional cell, index of negative top-dimensional cell)) + def cofaces_of_persistence_pairs(self): + """A persistence interval is described by a pair of cells, one that creates the + feature and one that kills it. The filtration values of those 2 cells give coordinates + for a point in a persistence diagram, or a bar in a barcode. Structurally, in the + cubical complexes provided here, the filtration value of any cell is the minimum of the + filtration values of the maximal cells that contain it. Connecting persistence diagram + coordinates to the corresponding value in the input (i.e. the filtration values of + the top-dimensional cells) is useful for differentiation purposes. + + This function returns a list of pairs of top-dimensional cells corresponding to + the persistence birth and death cells of the filtration. The cells are represented by + their indices in the input list of top-dimensional cells (and not their indices in the + internal datastructure that includes non-maximal cells). Note that when two adjacent + top-dimensional cells have the same filtration value, we arbitrarily return one of the two + when calling the function on one of their common faces. + + :returns: The top-dimensional cells/cofaces of the positive and negative cells. + :rtype: list of pairs(index of positive top-dimensional cell, index of negative top-dimensional cell) """ + cdef vector[vector[int]] persistence_result if self.pcohptr != NULL: - del self.pcohptr - self.pcohptr = new Cubical_complex_persistence_interface(self.thisptr, True) - cdef vector[pair[int, pair[pair[double, int], pair[double, int]]]] persistence_result - if self.pcohptr != NULL: - persistence_result = self.pcohptr.get_cofaces_of_cubical_persistence_pairs(homology_coeff_field, min_persistence) - return persistence_result + persistence_result = self.pcohptr.cofaces_of_cubical_persistence_pairs() + else: + print("cofaces_of_persistence_pairs function requires persistence function" + " to be launched first.") + return np.array(persistence_result) def betti_numbers(self): """This function returns the Betti numbers of the complex. diff --git a/src/python/include/Persistent_cohomology_interface.h b/src/python/include/Persistent_cohomology_interface.h index 1a1e716e..e5accf50 100644 --- a/src/python/include/Persistent_cohomology_interface.h +++ b/src/python/include/Persistent_cohomology_interface.h @@ -63,7 +63,6 @@ persistent_cohomology::Persistent_cohomology::get_persistent_pairs(); std::sort(std::begin(persistent_pairs), std::end(persistent_pairs), cmp); - std::vector>> persistence; for (auto pair : persistent_pairs) { persistence.push_back(std::make_pair(stptr_->dimension(get<0>(pair)), @@ -73,58 +72,52 @@ persistent_cohomology::Persistent_cohomology & cofaces, int splx){ - if (stptr_->dimension(stptr_->simplex(splx)) == stptr_->dimension()){cofaces.push_back(stptr_->simplex(splx));} - else{ for (auto v : stptr_->coboundary_simplex_range(stptr_->simplex(splx))){top_dimensional_cofaces(cofaces, stptr_->key(v));} } + int top_dimensional_coface(int splx){ + if (stptr_->dimension(splx) == stptr_->dimension()){return splx;} + else{ + for (auto v : stptr_->coboundary_simplex_range(splx)){ + if(stptr_->filtration(v) == stptr_->filtration(splx)){ + return top_dimensional_coface(v); + } + } + } } - std::vector, std::pair>>> get_cofaces_of_cubical_persistence_pairs(int homology_coeff_field, - double min_persistence) { + std::vector> cofaces_of_cubical_persistence_pairs() { // Warning: this function is meant to be used with CubicalComplex only!! + auto pairs = persistent_cohomology::Persistent_cohomology::get_persistent_pairs(); + // Gather all top-dimensional cells and store their simplex handles - std::vector max_splx; for (auto splx : stptr_->filtration_simplex_range()){ if (stptr_->dimension(splx) == stptr_->dimension()) max_splx.push_back(splx); } + std::vector max_splx; for (auto splx : stptr_->top_dimensional_cells_range()){ + max_splx.push_back(splx); + } // Sort these simplex handles and compute the ordering function // This function allows to go directly from the simplex handle to the position of the corresponding top-dimensional cell in the input data - std::map order; std::sort(max_splx.begin(), max_splx.end()); for (unsigned int i = 0; i < max_splx.size(); i++) order.insert(std::make_pair(max_splx[i], i)); - - persistent_cohomology::Persistent_cohomology::init_coefficients(homology_coeff_field); - persistent_cohomology::Persistent_cohomology::compute_persistent_cohomology(min_persistence); + std::map order; std::sort(max_splx.begin(), max_splx.end()); + for (unsigned int i = 0; i < max_splx.size(); i++) order.insert(std::make_pair(max_splx[i], i)); - // Custom sort and output persistence - cmp_intervals_by_dim_then_length cmp(stptr_); - auto persistent_pairs = persistent_cohomology::Persistent_cohomology::get_persistent_pairs(); - std::sort(std::begin(persistent_pairs), std::end(persistent_pairs), cmp); - - std::vector, std::pair>>> persistence; - for (auto pair : persistent_pairs) { - - double f0 = stptr_->filtration(get<0>(pair)); - // Recursively get the top-dimensional cells / cofaces associated to the persistence generator - std::vector faces0; top_dimensional_cofaces(faces0, stptr_->key(get<0>(pair))); - // Find the top-dimensional cell / coface with the same filtration value - int cf; for (unsigned int i = 0; i < faces0.size(); i++){if (stptr_->filtration(faces0[i]) == f0){cf = i; break;}} + std::vector> persistence_pairs; + for (auto pair : pairs) { + int h = stptr_->dimension(get<0>(pair)); + // Recursively get the top-dimensional cell / coface associated to the persistence generator + int face0 = top_dimensional_coface(get<0>(pair)); // Retrieve the index of the corresponding top-dimensional cell in the input data - int splx0 = order[faces0[cf]]; + int splx0 = order[face0]; int splx1 = -1; if (isfinite(stptr_->filtration(get<1>(pair)))){ - double f1 = stptr_->filtration(get<1>(pair)); - // Recursively get the top-dimensional cells / cofaces associated to the persistence generator - std::vector faces1; top_dimensional_cofaces(faces1, stptr_->key(get<1>(pair))); - // Find the top-dimensional cell / coface with the same filtration value - int cf; for (unsigned int i = 0; i < faces1.size(); i++){if (stptr_->filtration(faces1[i]) == f1){cf = i; break;}} + // Recursively get the top-dimensional cell / coface associated to the persistence generator + int face1 = top_dimensional_coface(get<1>(pair)); // Retrieve the index of the corresponding top-dimensional cell in the input data - splx1 = order[faces1[cf]]; + splx1 = order[face1]; } - - persistence.push_back(std::make_pair(stptr_->dimension(get<0>(pair)), std::make_pair(std::make_pair(stptr_->filtration(get<0>(pair)), splx0), std::make_pair(stptr_->filtration(get<1>(pair)), splx1)))); + std::vector vect{ h, splx0, splx1}; + persistence_pairs.push_back(vect); } - return persistence; + return persistence_pairs; } std::vector, std::vector>> persistence_pairs() { -- cgit v1.2.3 From 5a737eefc7abd690e8a174d2557d0157e77f5f4c Mon Sep 17 00:00:00 2001 From: mathieu Date: Tue, 10 Mar 2020 19:13:37 -0400 Subject: new fixes --- .../include/gudhi/Bitmap_cubical_complex.h | 1 - src/python/gudhi/cubical_complex.pyx | 32 +++++++++++----------- src/python/gudhi/periodic_cubical_complex.pyx | 28 +++++++++++++++++++ src/python/gudhi/simplex_tree.pyx | 2 +- .../include/Persistent_cohomology_interface.h | 3 +- src/python/test/test_cubical_complex.py | 5 ++++ 6 files changed, 52 insertions(+), 19 deletions(-) (limited to 'src') diff --git a/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex.h b/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex.h index bf09532e..37514dee 100644 --- a/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex.h +++ b/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex.h @@ -340,7 +340,6 @@ class Bitmap_cubical_complex : public T { * that provides ranges for the Boundary_simplex_iterator. **/ Boundary_simplex_range boundary_simplex_range(Simplex_handle sh) { return this->get_boundary_of_a_cell(sh); } - Boundary_simplex_range coboundary_simplex_range(Simplex_handle sh) { return this->get_coboundary_of_a_cell(sh); } /** * filtration_simplex_range creates an object of a Filtration_simplex_range class diff --git a/src/python/gudhi/cubical_complex.pyx b/src/python/gudhi/cubical_complex.pyx index 8cf43539..9e701fe6 100644 --- a/src/python/gudhi/cubical_complex.pyx +++ b/src/python/gudhi/cubical_complex.pyx @@ -148,22 +148,22 @@ cdef class CubicalComplex: def cofaces_of_persistence_pairs(self): """A persistence interval is described by a pair of cells, one that creates the - feature and one that kills it. The filtration values of those 2 cells give coordinates - for a point in a persistence diagram, or a bar in a barcode. Structurally, in the - cubical complexes provided here, the filtration value of any cell is the minimum of the - filtration values of the maximal cells that contain it. Connecting persistence diagram - coordinates to the corresponding value in the input (i.e. the filtration values of - the top-dimensional cells) is useful for differentiation purposes. - - This function returns a list of pairs of top-dimensional cells corresponding to - the persistence birth and death cells of the filtration. The cells are represented by - their indices in the input list of top-dimensional cells (and not their indices in the - internal datastructure that includes non-maximal cells). Note that when two adjacent - top-dimensional cells have the same filtration value, we arbitrarily return one of the two - when calling the function on one of their common faces. - - :returns: The top-dimensional cells/cofaces of the positive and negative cells. - :rtype: list of pairs(index of positive top-dimensional cell, index of negative top-dimensional cell) + feature and one that kills it. The filtration values of those 2 cells give coordinates + for a point in a persistence diagram, or a bar in a barcode. Structurally, in the + cubical complexes provided here, the filtration value of any cell is the minimum of the + filtration values of the maximal cells that contain it. Connecting persistence diagram + coordinates to the corresponding value in the input (i.e. the filtration values of + the top-dimensional cells) is useful for differentiation purposes. + + This function returns a list of pairs of top-dimensional cells corresponding to + the persistence birth and death cells of the filtration. The cells are represented by + their indices in the input list of top-dimensional cells (and not their indices in the + internal datastructure that includes non-maximal cells). Note that when two adjacent + top-dimensional cells have the same filtration value, we arbitrarily return one of the two + when calling the function on one of their common faces. + + :returns: The top-dimensional cells/cofaces of the positive and negative cells, together with the corresponding homological dimension. + :rtype: numpy array of integers of shape [number_of_persistence_points, 3], the integers of eah row being: (homological dimension, index of positive top-dimensional cell, index of negative top-dimensional cell). If the homological feature is essential, i.e., if the death time is +infinity, then the index of the corresponding negative top-dimensional cell is -1. """ cdef vector[vector[int]] persistence_result if self.pcohptr != NULL: diff --git a/src/python/gudhi/periodic_cubical_complex.pyx b/src/python/gudhi/periodic_cubical_complex.pyx index 37f76201..ba039e80 100644 --- a/src/python/gudhi/periodic_cubical_complex.pyx +++ b/src/python/gudhi/periodic_cubical_complex.pyx @@ -31,6 +31,7 @@ cdef extern from "Persistent_cohomology_interface.h" namespace "Gudhi": cdef cppclass Periodic_cubical_complex_persistence_interface "Gudhi::Persistent_cohomology_interface>>": Periodic_cubical_complex_persistence_interface(Periodic_cubical_complex_base_interface * st, bool persistence_dim_max) vector[pair[int, pair[double, double]]] get_persistence(int homology_coeff_field, double min_persistence) + vector[vector[int]] cofaces_of_cubical_persistence_pairs() vector[int] betti_numbers() vector[int] persistent_betti_numbers(double from_value, double to_value) vector[pair[double,double]] intervals_in_dimension(int dimension) @@ -155,6 +156,33 @@ cdef class PeriodicCubicalComplex: persistence_result = self.pcohptr.get_persistence(homology_coeff_field, min_persistence) return persistence_result + def cofaces_of_persistence_pairs(self): + """A persistence interval is described by a pair of cells, one that creates the + feature and one that kills it. The filtration values of those 2 cells give coordinates + for a point in a persistence diagram, or a bar in a barcode. Structurally, in the + cubical complexes provided here, the filtration value of any cell is the minimum of the + filtration values of the maximal cells that contain it. Connecting persistence diagram + coordinates to the corresponding value in the input (i.e. the filtration values of + the top-dimensional cells) is useful for differentiation purposes. + + This function returns a list of pairs of top-dimensional cells corresponding to + the persistence birth and death cells of the filtration. The cells are represented by + their indices in the input list of top-dimensional cells (and not their indices in the + internal datastructure that includes non-maximal cells). Note that when two adjacent + top-dimensional cells have the same filtration value, we arbitrarily return one of the two + when calling the function on one of their common faces. + + :returns: The top-dimensional cells/cofaces of the positive and negative cells, together with the corresponding homological dimension. + :rtype: numpy array of integers of shape [number_of_persistence_points, 3], the integers of eah row being: (homological dimension, index of positive top-dimensional cell, index of negative top-dimensional cell). If the homological feature is essential, i.e., if the death time is +infinity, then the index of the corresponding negative top-dimensional cell is -1. + """ + cdef vector[vector[int]] persistence_result + if self.pcohptr != NULL: + persistence_result = self.pcohptr.cofaces_of_cubical_persistence_pairs() + else: + print("cofaces_of_persistence_pairs function requires persistence function" + " to be launched first.") + return np.array(persistence_result) + def betti_numbers(self): """This function returns the Betti numbers of the complex. diff --git a/src/python/gudhi/simplex_tree.pyx b/src/python/gudhi/simplex_tree.pyx index 85d25492..b18627c4 100644 --- a/src/python/gudhi/simplex_tree.pyx +++ b/src/python/gudhi/simplex_tree.pyx @@ -508,7 +508,7 @@ cdef class SimplexTree: """ if self.pcohptr != NULL: if persistence_file != '': - self.pcohptr.write_output_diagram(str.encode(persistence_file)) + self.pcohptr.write_output_diagram(persistence_file.encode('utf-8')) else: print("persistence_file must be specified") else: diff --git a/src/python/include/Persistent_cohomology_interface.h b/src/python/include/Persistent_cohomology_interface.h index e5accf50..defac88c 100644 --- a/src/python/include/Persistent_cohomology_interface.h +++ b/src/python/include/Persistent_cohomology_interface.h @@ -75,12 +75,13 @@ persistent_cohomology::Persistent_cohomologydimension(splx) == stptr_->dimension()){return splx;} else{ - for (auto v : stptr_->coboundary_simplex_range(splx)){ + for (auto v : stptr_->get_coboundary_of_a_cell(splx)){ if(stptr_->filtration(v) == stptr_->filtration(splx)){ return top_dimensional_coface(v); } } } + return splx; } std::vector> cofaces_of_cubical_persistence_pairs() { diff --git a/src/python/test/test_cubical_complex.py b/src/python/test/test_cubical_complex.py index 8c1b2600..8af63355 100755 --- a/src/python/test/test_cubical_complex.py +++ b/src/python/test/test_cubical_complex.py @@ -147,3 +147,8 @@ def test_connected_sublevel_sets(): periodic_dimensions = periodic_dimensions) assert cub.persistence() == [(0, (2.0, float("inf")))] assert cub.betti_numbers() == [1, 0, 0] + +def test_connected_sublevel_sets(): + cub = CubicalComplex(top_dimensional_cells = [[0, 0, 0], [0, 1, 0], [0, 0, 0]]) + cub.persistence() + assert cub.cofaces_of_persistence_pairs() == np.array([[1, 7, 4], [0, 8, -1]]) -- cgit v1.2.3 From 45b918a17cfa26a0c58d7871b869aa13b0e45019 Mon Sep 17 00:00:00 2001 From: mathieu Date: Wed, 11 Mar 2020 12:05:15 -0400 Subject: moved location of top_dimensional_coface function --- ext/hera | 2 +- .../include/gudhi/Bitmap_cubical_complex_base.h | 21 +++++++++++++++++++++ src/python/gudhi/cubical_complex.pyx | 4 +++- src/python/gudhi/periodic_cubical_complex.pyx | 4 +++- .../include/Persistent_cohomology_interface.h | 16 ++-------------- 5 files changed, 30 insertions(+), 17 deletions(-) (limited to 'src') diff --git a/ext/hera b/ext/hera index cb1838e6..9a899718 160000 --- a/ext/hera +++ b/ext/hera @@ -1 +1 @@ -Subproject commit cb1838e682ec07f80720241cf9098400caeb83c7 +Subproject commit 9a89971855acefe39dce0e2adadf53b88ca8f683 diff --git a/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_base.h b/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_base.h index 0d6299d2..7496d74a 100644 --- a/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_base.h +++ b/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_base.h @@ -109,6 +109,14 @@ class Bitmap_cubical_complex_base { **/ virtual inline std::vector get_coboundary_of_a_cell(std::size_t cell) const; + /** + * This function computes the index of one of the top-dimensional cubes (chosen arbitrarily) associated + * to a given simplex handle. Note that the input parameter is not necessarily a cube, it might also + * be an edge or vertex of a cube. On the other hand, the output is always indicating the position of + * a cube in the data structure. + **/ + inline int get_top_dimensional_coface_of_a_cell(int splx); + /** * This procedure compute incidence numbers between cubes. For a cube \f$A\f$ of * dimension n and a cube \f$B \subset A\f$ of dimension n-1, an incidence @@ -602,6 +610,19 @@ void Bitmap_cubical_complex_base::setup_bitmap_based_on_top_dimensional_cells this->impose_lower_star_filtration(); } +template +int Bitmap_cubical_complex_base::get_top_dimensional_coface_of_a_cell(int splx) { + if (this->get_dimension_of_a_cell(splx) == this->dimension()){return splx;} + else{ + for (auto v : this->get_coboundary_of_a_cell(splx)){ + if(this->get_cell_data(v) == this->get_cell_data(splx)){ + return this->get_top_dimensional_coface_of_a_cell(v); + } + } + } + return splx; +} + template Bitmap_cubical_complex_base::Bitmap_cubical_complex_base(const std::vector& sizes_in_following_directions, const std::vector& top_dimensional_cells) { diff --git a/src/python/gudhi/cubical_complex.pyx b/src/python/gudhi/cubical_complex.pyx index 9e701fe6..84fec60e 100644 --- a/src/python/gudhi/cubical_complex.pyx +++ b/src/python/gudhi/cubical_complex.pyx @@ -163,7 +163,9 @@ cdef class CubicalComplex: when calling the function on one of their common faces. :returns: The top-dimensional cells/cofaces of the positive and negative cells, together with the corresponding homological dimension. - :rtype: numpy array of integers of shape [number_of_persistence_points, 3], the integers of eah row being: (homological dimension, index of positive top-dimensional cell, index of negative top-dimensional cell). If the homological feature is essential, i.e., if the death time is +infinity, then the index of the corresponding negative top-dimensional cell is -1. + :rtype: numpy array of integers of shape [number_of_persistence_points, 3], the integers of eah row being: (homological dimension, + index of positive top-dimensional cell, index of negative top-dimensional cell). If the homological feature is essential, i.e., if + the death time is +infinity, then the index of the corresponding negative top-dimensional cell is -1. """ cdef vector[vector[int]] persistence_result if self.pcohptr != NULL: diff --git a/src/python/gudhi/periodic_cubical_complex.pyx b/src/python/gudhi/periodic_cubical_complex.pyx index ba039e80..993d95c7 100644 --- a/src/python/gudhi/periodic_cubical_complex.pyx +++ b/src/python/gudhi/periodic_cubical_complex.pyx @@ -173,7 +173,9 @@ cdef class PeriodicCubicalComplex: when calling the function on one of their common faces. :returns: The top-dimensional cells/cofaces of the positive and negative cells, together with the corresponding homological dimension. - :rtype: numpy array of integers of shape [number_of_persistence_points, 3], the integers of eah row being: (homological dimension, index of positive top-dimensional cell, index of negative top-dimensional cell). If the homological feature is essential, i.e., if the death time is +infinity, then the index of the corresponding negative top-dimensional cell is -1. + :rtype: numpy array of integers of shape [number_of_persistence_points, 3], the integers of eah row being: (homological dimension, + index of positive top-dimensional cell, index of negative top-dimensional cell). If the homological feature is essential, i.e., if + the death time is +infinity, then the index of the corresponding negative top-dimensional cell is -1. """ cdef vector[vector[int]] persistence_result if self.pcohptr != NULL: diff --git a/src/python/include/Persistent_cohomology_interface.h b/src/python/include/Persistent_cohomology_interface.h index defac88c..77555349 100644 --- a/src/python/include/Persistent_cohomology_interface.h +++ b/src/python/include/Persistent_cohomology_interface.h @@ -72,18 +72,6 @@ persistent_cohomology::Persistent_cohomologydimension(splx) == stptr_->dimension()){return splx;} - else{ - for (auto v : stptr_->get_coboundary_of_a_cell(splx)){ - if(stptr_->filtration(v) == stptr_->filtration(splx)){ - return top_dimensional_coface(v); - } - } - } - return splx; - } - std::vector> cofaces_of_cubical_persistence_pairs() { // Warning: this function is meant to be used with CubicalComplex only!! @@ -104,14 +92,14 @@ persistent_cohomology::Persistent_cohomologydimension(get<0>(pair)); // Recursively get the top-dimensional cell / coface associated to the persistence generator - int face0 = top_dimensional_coface(get<0>(pair)); + int face0 = stptr_->get_top_dimensional_coface_of_a_cell(get<0>(pair)); // Retrieve the index of the corresponding top-dimensional cell in the input data int splx0 = order[face0]; int splx1 = -1; if (isfinite(stptr_->filtration(get<1>(pair)))){ // Recursively get the top-dimensional cell / coface associated to the persistence generator - int face1 = top_dimensional_coface(get<1>(pair)); + int face1 = stptr_->get_top_dimensional_coface_of_a_cell(get<1>(pair)); // Retrieve the index of the corresponding top-dimensional cell in the input data splx1 = order[face1]; } -- cgit v1.2.3 From 5ce1ee8976ced78de839ef629522c95324b2fabd Mon Sep 17 00:00:00 2001 From: yuichi-ike Date: Mon, 6 Apr 2020 16:25:27 +0900 Subject: weighted rips added --- src/python/CMakeLists.txt | 3 +++ src/python/gudhi/weighted_rips_complex.py | 41 +++++++++++++++++++++++++++++++ 2 files changed, 44 insertions(+) create mode 100644 src/python/gudhi/weighted_rips_complex.py (limited to 'src') diff --git a/src/python/CMakeLists.txt b/src/python/CMakeLists.txt index d7a6a4db..cac4553a 100644 --- a/src/python/CMakeLists.txt +++ b/src/python/CMakeLists.txt @@ -415,6 +415,9 @@ if(PYTHONINTERP_FOUND) add_gudhi_py_test(test_dtm) endif() + # Weighted Rips + add_gudhi_py_test(test_weighted_rips) + # Documentation generation is available through sphinx - requires all modules if(SPHINX_PATH) if(MATPLOTLIB_FOUND) diff --git a/src/python/gudhi/weighted_rips_complex.py b/src/python/gudhi/weighted_rips_complex.py new file mode 100644 index 00000000..34a627cb --- /dev/null +++ b/src/python/gudhi/weighted_rips_complex.py @@ -0,0 +1,41 @@ +# 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): Raphaël Tinarrage and Yuichi Ike +# +# Copyright (C) 2020 Inria, Copyright (C) 2020 FUjitsu Laboratories Ltd. +# +# Modification(s): +# - YYYY/MM Author: Description of the modification + +from gudhi import SimplexTree + +class WeightedRipsComplex: + """ + class to generate a weighted Rips complex + from a distance matrix and filtration value + """ + def __init__(self, + distance_matrix=None, + filtration_values=None, + max_filtration=float('inf'), sparse=None): + self.distance_matrix = distance_matrix + self.filtration_values = filtration_values + self.max_filtration = max_filtration + + def create_simplex_tree(self, max_dimension): + dist = self.distance_matrix + F = self.filtration_values + num_pts = len(dist) + + st = SimplexTree() + + for i in range(num_pts): + if F[i] < self.max_filtration: + st.insert([i], F[i]) + for i in range(num_pts): + for j in range(num_pts): + value = (dist[i][j] + F[i] + F[j]) / 2 + if value < self.max_filtration: + st.insert([i,j], filtration=value) + return st + -- cgit v1.2.3 From fadeb80b46001779e2a998941a02195921b03124 Mon Sep 17 00:00:00 2001 From: yuichi-ike Date: Mon, 6 Apr 2020 16:31:59 +0900 Subject: test_weighted_rips added --- src/python/test/test_weighted_rips.py | 27 +++++++++++++++++++++++++++ 1 file changed, 27 insertions(+) create mode 100644 src/python/test/test_weighted_rips.py (limited to 'src') diff --git a/src/python/test/test_weighted_rips.py b/src/python/test/test_weighted_rips.py new file mode 100644 index 00000000..f0db6798 --- /dev/null +++ b/src/python/test/test_weighted_rips.py @@ -0,0 +1,27 @@ +""" 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): Yuichi Ike + + Copyright (C) 2020 Inria + + Modification(s): + - YYYY/MM Author: Description of the modification +""" + +from gudhi.weighted_rips_complex import WeightedRipsComplex +from gudhi.point_cloud.dtm import DTM +import numpy +from scipy.spatial.distance import cdist +import pytest + +def test_dtm_rips_complex(): + pts = numpy.array([[2.0, 2], [0, 1], [3, 4]]) + dist = cdist(pts,pts) + dtm = DTM(2, q=2, metric="precomputed") + r = dtm.fit_transform(dist) + w_rips = WeightedRipsComplex(distance_mattix=dist, filtration_values=r) + st = w_rips.create_simplex_tree(max_dimension=2) + diag = st.persistence() + assert diag == [(0, (1.5811388300841898, float("inf"))), (0, (1.5811388300841898, 2.699172818834085)), (0, (1.5811388300841898, 2.699172818834085))] + + \ No newline at end of file -- cgit v1.2.3 From 5737c5e1e89cc4c939a784742f25b26ca163332d Mon Sep 17 00:00:00 2001 From: yuichi-ike Date: Mon, 6 Apr 2020 16:43:55 +0900 Subject: comments added --- src/python/gudhi/weighted_rips_complex.py | 16 +++++++++++++++- 1 file changed, 15 insertions(+), 1 deletion(-) (limited to 'src') diff --git a/src/python/gudhi/weighted_rips_complex.py b/src/python/gudhi/weighted_rips_complex.py index 34a627cb..84e8e38e 100644 --- a/src/python/gudhi/weighted_rips_complex.py +++ b/src/python/gudhi/weighted_rips_complex.py @@ -17,12 +17,26 @@ class WeightedRipsComplex: def __init__(self, distance_matrix=None, filtration_values=None, - max_filtration=float('inf'), sparse=None): + max_filtration=float('inf')): + """ + Parameters: + distance_matrix: list of list of float, + distance matrix (full square or lower triangular) + filtration_values: list of float, + flitration value for each index + max_filtration: float, + specifies the maximal filtration value to be considered + """ self.distance_matrix = distance_matrix self.filtration_values = filtration_values self.max_filtration = max_filtration def create_simplex_tree(self, max_dimension): + """ + Parameter: + max_dimension: int + graph expansion until this given dimension + """ dist = self.distance_matrix F = self.filtration_values num_pts = len(dist) -- cgit v1.2.3 From 15586d479be885319dde6f703c3126176b796732 Mon Sep 17 00:00:00 2001 From: yuichi-ike Date: Mon, 6 Apr 2020 16:48:21 +0900 Subject: bug fixed --- src/python/gudhi/weighted_rips_complex.py | 2 ++ 1 file changed, 2 insertions(+) (limited to 'src') diff --git a/src/python/gudhi/weighted_rips_complex.py b/src/python/gudhi/weighted_rips_complex.py index 84e8e38e..7d14ac65 100644 --- a/src/python/gudhi/weighted_rips_complex.py +++ b/src/python/gudhi/weighted_rips_complex.py @@ -51,5 +51,7 @@ class WeightedRipsComplex: value = (dist[i][j] + F[i] + F[j]) / 2 if value < self.max_filtration: st.insert([i,j], filtration=value) + + st.expansion(max_dimension) return st -- cgit v1.2.3 From a4fa5f673784a842e9fac13003c843d454c888a4 Mon Sep 17 00:00:00 2001 From: yuichi-ike Date: Mon, 6 Apr 2020 21:19:55 +0900 Subject: bug fixed, parameter name changed --- src/python/CMakeLists.txt | 2 ++ src/python/gudhi/weighted_rips_complex.py | 19 +++++++++++-------- src/python/test/test_weighted_rips.py | 13 ++++++------- 3 files changed, 19 insertions(+), 15 deletions(-) (limited to 'src') diff --git a/src/python/CMakeLists.txt b/src/python/CMakeLists.txt index cac4553a..4b87ed9b 100644 --- a/src/python/CMakeLists.txt +++ b/src/python/CMakeLists.txt @@ -57,6 +57,7 @@ if(PYTHONINTERP_FOUND) set(GUDHI_PYTHON_MODULES_EXTRA "${GUDHI_PYTHON_MODULES_EXTRA}'representations', ") set(GUDHI_PYTHON_MODULES_EXTRA "${GUDHI_PYTHON_MODULES_EXTRA}'wasserstein', ") set(GUDHI_PYTHON_MODULES_EXTRA "${GUDHI_PYTHON_MODULES_EXTRA}'point_cloud', ") + set(GUDHI_PYTHON_MODULES_EXTRA "${GUDHI_PYTHON_MODULES_EXTRA}'weighted_rips_complex', ") add_gudhi_debug_info("Python version ${PYTHON_VERSION_STRING}") add_gudhi_debug_info("Cython version ${CYTHON_VERSION}") @@ -228,6 +229,7 @@ if(PYTHONINTERP_FOUND) file(COPY "gudhi/representations" DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/gudhi/") file(COPY "gudhi/wasserstein.py" DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/gudhi") file(COPY "gudhi/point_cloud" DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/gudhi") + file(COPY "gudhi/weighted_rips_complex.py" DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/gudhi") add_custom_command( OUTPUT gudhi.so diff --git a/src/python/gudhi/weighted_rips_complex.py b/src/python/gudhi/weighted_rips_complex.py index 7d14ac65..9df2ddf9 100644 --- a/src/python/gudhi/weighted_rips_complex.py +++ b/src/python/gudhi/weighted_rips_complex.py @@ -1,6 +1,6 @@ # 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): Raphaël Tinarrage and Yuichi Ike +# Author(s): Raphaël Tinarrage, Yuichi Ike, Masatoshi Takenouchi # # Copyright (C) 2020 Inria, Copyright (C) 2020 FUjitsu Laboratories Ltd. # @@ -12,23 +12,26 @@ from gudhi import SimplexTree class WeightedRipsComplex: """ class to generate a weighted Rips complex - from a distance matrix and filtration value + from a distance matrix and weights on vertices """ def __init__(self, - distance_matrix=None, - filtration_values=None, + distance_matrix, + weights=None, max_filtration=float('inf')): """ Parameters: distance_matrix: list of list of float, distance matrix (full square or lower triangular) filtration_values: list of float, - flitration value for each index + weight for each vertex max_filtration: float, specifies the maximal filtration value to be considered """ self.distance_matrix = distance_matrix - self.filtration_values = filtration_values + if weights is not None: + self.weights = weights + else: + self.weights = [0] * len(distance_matrix) self.max_filtration = max_filtration def create_simplex_tree(self, max_dimension): @@ -38,7 +41,7 @@ class WeightedRipsComplex: graph expansion until this given dimension """ dist = self.distance_matrix - F = self.filtration_values + F = self.weights num_pts = len(dist) st = SimplexTree() @@ -47,7 +50,7 @@ class WeightedRipsComplex: if F[i] < self.max_filtration: st.insert([i], F[i]) for i in range(num_pts): - for j in range(num_pts): + for j in range(i): value = (dist[i][j] + F[i] + F[j]) / 2 if value < self.max_filtration: st.insert([i,j], filtration=value) diff --git a/src/python/test/test_weighted_rips.py b/src/python/test/test_weighted_rips.py index f0db6798..7896fb78 100644 --- a/src/python/test/test_weighted_rips.py +++ b/src/python/test/test_weighted_rips.py @@ -1,6 +1,6 @@ """ 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): Yuichi Ike + Author(s): Yuichi Ike and Masatoshi Takenouchi Copyright (C) 2020 Inria @@ -10,18 +10,17 @@ from gudhi.weighted_rips_complex import WeightedRipsComplex from gudhi.point_cloud.dtm import DTM -import numpy +import numpy as np from scipy.spatial.distance import cdist import pytest def test_dtm_rips_complex(): - pts = numpy.array([[2.0, 2], [0, 1], [3, 4]]) + pts = np.array([[2.0, 2], [0, 1], [3, 4]]) dist = cdist(pts,pts) dtm = DTM(2, q=2, metric="precomputed") r = dtm.fit_transform(dist) - w_rips = WeightedRipsComplex(distance_mattix=dist, filtration_values=r) + w_rips = WeightedRipsComplex(distance_mattix=dist, weights=r) st = w_rips.create_simplex_tree(max_dimension=2) - diag = st.persistence() - assert diag == [(0, (1.5811388300841898, float("inf"))), (0, (1.5811388300841898, 2.699172818834085)), (0, (1.5811388300841898, 2.699172818834085))] + persistence_intervals0 = st.persistence_intervals_in_dimension(0) + assert persistence_intervals0 == pytest.approx(np.array([[1.58113883, 2.69917282],[1.58113883, 2.69917282], [1.58113883, float("inf")]])) - \ No newline at end of file -- cgit v1.2.3 From 4294e5fc6e1bff246a7d22f1bd98f91b62f14163 Mon Sep 17 00:00:00 2001 From: yuichi-ike Date: Tue, 7 Apr 2020 09:36:03 +0900 Subject: filtration value fixed --- src/python/gudhi/weighted_rips_complex.py | 2 +- src/python/test/test_weighted_rips.py | 12 +++++++++++- 2 files changed, 12 insertions(+), 2 deletions(-) (limited to 'src') diff --git a/src/python/gudhi/weighted_rips_complex.py b/src/python/gudhi/weighted_rips_complex.py index 9df2ddf9..7e504b2c 100644 --- a/src/python/gudhi/weighted_rips_complex.py +++ b/src/python/gudhi/weighted_rips_complex.py @@ -51,7 +51,7 @@ class WeightedRipsComplex: st.insert([i], F[i]) for i in range(num_pts): for j in range(i): - value = (dist[i][j] + F[i] + F[j]) / 2 + value = max(F[i], F[j], (dist[i][j] + F[i] + F[j]) / 2) if value < self.max_filtration: st.insert([i,j], filtration=value) diff --git a/src/python/test/test_weighted_rips.py b/src/python/test/test_weighted_rips.py index 7896fb78..a3235276 100644 --- a/src/python/test/test_weighted_rips.py +++ b/src/python/test/test_weighted_rips.py @@ -14,13 +14,23 @@ import numpy as np from scipy.spatial.distance import cdist import pytest +def test_non_dtm_rips_complex(): + dist = [[], [1]] + weights = [1, 100] + w_rips = WeightedRipsComplex(distance_matrix=dist, weights=weights) + st = w_rips.create_simplex_tree(max_dimension=2) + assert st.filtration([0,1]) == pytest.approx(100.0) + + def test_dtm_rips_complex(): pts = np.array([[2.0, 2], [0, 1], [3, 4]]) dist = cdist(pts,pts) dtm = DTM(2, q=2, metric="precomputed") r = dtm.fit_transform(dist) - w_rips = WeightedRipsComplex(distance_mattix=dist, weights=r) + w_rips = WeightedRipsComplex(distance_matrix=dist, weights=r) st = w_rips.create_simplex_tree(max_dimension=2) + st.persistence() persistence_intervals0 = st.persistence_intervals_in_dimension(0) assert persistence_intervals0 == pytest.approx(np.array([[1.58113883, 2.69917282],[1.58113883, 2.69917282], [1.58113883, float("inf")]])) + -- cgit v1.2.3 From c2b6d95f0b01ca913ddc704350cbfe37bcf13c3a Mon Sep 17 00:00:00 2001 From: MathieuCarriere Date: Tue, 28 Apr 2020 19:28:24 -0400 Subject: update output --- .../include/gudhi/Bitmap_cubical_complex_base.h | 5 ++-- src/python/gudhi/cubical_complex.pyx | 33 ++++++++++++++++++---- src/python/gudhi/periodic_cubical_complex.pyx | 33 ++++++++++++++++++---- .../include/Persistent_cohomology_interface.h | 6 ++-- src/python/test/test_cubical_complex.py | 7 ++++- 5 files changed, 69 insertions(+), 15 deletions(-) (limited to 'src') diff --git a/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_base.h b/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_base.h index 6441c129..248ebdb6 100644 --- a/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_base.h +++ b/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_base.h @@ -110,8 +110,9 @@ class Bitmap_cubical_complex_base { virtual inline std::vector get_coboundary_of_a_cell(std::size_t cell) const; /** - * This function computes the index of one of the top-dimensional cubes (chosen arbitrarily) associated - * to a given simplex handle. Note that the input parameter is not necessarily a cube, it might also + * This function finds a top-dimensional cell that is incident to the input cell and has + * the same filtration value. In case several cells are suitable, an arbitrary one is + * returned. Note that the input parameter is not necessarily a cube, it might also * be an edge or vertex of a cube. On the other hand, the output is always indicating the position of * a cube in the data structure. **/ diff --git a/src/python/gudhi/cubical_complex.pyx b/src/python/gudhi/cubical_complex.pyx index 69d0f0b6..884b0664 100644 --- a/src/python/gudhi/cubical_complex.pyx +++ b/src/python/gudhi/cubical_complex.pyx @@ -187,18 +187,41 @@ cdef class CubicalComplex: top-dimensional cells have the same filtration value, we arbitrarily return one of the two when calling the function on one of their common faces. - :returns: The top-dimensional cells/cofaces of the positive and negative cells, together with the corresponding homological dimension. - :rtype: numpy array of integers of shape [number_of_persistence_points, 3], the integers of eah row being: (homological dimension, - index of positive top-dimensional cell, index of negative top-dimensional cell). If the homological feature is essential, i.e., if - the death time is +infinity, then the index of the corresponding negative top-dimensional cell is -1. + :returns: The top-dimensional cells/cofaces of the positive and negative cells, + together with the corresponding homological dimension, in two lists of numpy arrays of integers. + The first list contains the regular persistence pairs, grouped by dimension. + It contains numpy arrays of shape [number_of_persistence_points, 2]. + The indices of the arrays in the list correspond to the homological dimensions, and the + integers of each row in each array correspond to: (index of positive top-dimensional cell, + index of negative top-dimensional cell). + The second list contains the essential features, grouped by dimension. + It contains numpy arrays of shape [number_of_persistence_points, 1]. + The indices of the arrays in the list correspond to the homological dimensions, and the + integers of each row in each array correspond to: (index of positive top-dimensional cell). """ cdef vector[vector[int]] persistence_result if self.pcohptr != NULL: + output = [[],[]] persistence_result = self.pcohptr.cofaces_of_cubical_persistence_pairs() + pr = np.array(persistence_result) + + ess_ind = np.argwhere(pr[:,2] == -1)[:,0] + ess = pr[ess_ind] + max_h = max(ess[:,0])+1 + for h in range(max_h): + hidxs = np.argwhere(ess[:,0] == h)[:,0] + output[1].append(ess[hidxs][:,1]) + + reg_ind = np.setdiff1d(np.array(range(len(pr))), ess_ind) + reg = pr[reg_ind] + max_h = max(reg[:,0])+1 + for h in range(max_h): + hidxs = np.argwhere(reg[:,0] == h)[:,0] + output[0].append(reg[hidxs][:,1:]) else: print("cofaces_of_persistence_pairs function requires persistence function" " to be launched first.") - return np.array(persistence_result) + return output def betti_numbers(self): """This function returns the Betti numbers of the complex. diff --git a/src/python/gudhi/periodic_cubical_complex.pyx b/src/python/gudhi/periodic_cubical_complex.pyx index 78565cf8..3cf2ff01 100644 --- a/src/python/gudhi/periodic_cubical_complex.pyx +++ b/src/python/gudhi/periodic_cubical_complex.pyx @@ -192,18 +192,41 @@ cdef class PeriodicCubicalComplex: top-dimensional cells have the same filtration value, we arbitrarily return one of the two when calling the function on one of their common faces. - :returns: The top-dimensional cells/cofaces of the positive and negative cells, together with the corresponding homological dimension. - :rtype: numpy array of integers of shape [number_of_persistence_points, 3], the integers of eah row being: (homological dimension, - index of positive top-dimensional cell, index of negative top-dimensional cell). If the homological feature is essential, i.e., if - the death time is +infinity, then the index of the corresponding negative top-dimensional cell is -1. + :returns: The top-dimensional cells/cofaces of the positive and negative cells, + together with the corresponding homological dimension, in two lists of numpy arrays of integers. + The first list contains the regular persistence pairs, grouped by dimension. + It contains numpy arrays of shape [number_of_persistence_points, 2]. + The indices of the arrays in the list correspond to the homological dimensions, and the + integers of each row in each array correspond to: (index of positive top-dimensional cell, + index of negative top-dimensional cell). + The second list contains the essential features, grouped by dimension. + It contains numpy arrays of shape [number_of_persistence_points, 1]. + The indices of the arrays in the list correspond to the homological dimensions, and the + integers of each row in each array correspond to: (index of positive top-dimensional cell). """ cdef vector[vector[int]] persistence_result if self.pcohptr != NULL: + output = [[],[]] persistence_result = self.pcohptr.cofaces_of_cubical_persistence_pairs() + pr = np.array(persistence_result) + + ess_ind = np.argwhere(pr[:,2] == -1)[:,0] + ess = pr[ess_ind] + max_h = max(ess[:,0])+1 + for h in range(max_h): + hidxs = np.argwhere(ess[:,0] == h)[:,0] + output[1].append(ess[hidxs][:,1]) + + reg_ind = np.setdiff1d(np.array(range(len(pr))), ess_ind) + reg = pr[reg_ind] + max_h = max(reg[:,0])+1 + for h in range(max_h): + hidxs = np.argwhere(reg[:,0] == h)[:,0] + output[0].append(reg[hidxs][:,1:]) else: print("cofaces_of_persistence_pairs function requires persistence function" " to be launched first.") - return np.array(persistence_result) + return output def betti_numbers(self): """This function returns the Betti numbers of the complex. diff --git a/src/python/include/Persistent_cohomology_interface.h b/src/python/include/Persistent_cohomology_interface.h index 59024212..32e6ee9c 100644 --- a/src/python/include/Persistent_cohomology_interface.h +++ b/src/python/include/Persistent_cohomology_interface.h @@ -16,6 +16,7 @@ #include #include // for std::pair #include // for sort +#include namespace Gudhi { @@ -80,8 +81,9 @@ persistent_cohomology::Persistent_cohomology order; std::sort(max_splx.begin(), max_splx.end()); - for (unsigned int i = 0; i < max_splx.size(); i++) order.insert(std::make_pair(max_splx[i], i)); + std::unordered_map order; + //std::sort(max_splx.begin(), max_splx.end()); + for (unsigned int i = 0; i < max_splx.size(); i++) order.emplace(max_splx[i], i); std::vector> persistence_pairs; for (auto pair : pairs) { diff --git a/src/python/test/test_cubical_complex.py b/src/python/test/test_cubical_complex.py index dd7653c2..5c59db8f 100755 --- a/src/python/test/test_cubical_complex.py +++ b/src/python/test/test_cubical_complex.py @@ -151,4 +151,9 @@ def test_connected_sublevel_sets(): def test_cubical_generators(): cub = CubicalComplex(top_dimensional_cells = [[0, 0, 0], [0, 1, 0], [0, 0, 0]]) cub.persistence() - assert np.array_equal(cub.cofaces_of_persistence_pairs(), np.array([[1, 7, 4], [0, 8, -1]])) + g = cub.cofaces_of_persistence_pairs() + assert len(g[0]) == 2 + assert len(g[1]) == 1 + assert np.array_equal(g[0][0], np.empty(shape=[0,2])) + assert np.array_equal(g[0][1], np.array([[7, 4]])) + assert np.array_equal(g[1][0], np.array([8])) -- cgit v1.2.3 From 2b5586fd60848b159fb4fa4481e61bab0e0cd766 Mon Sep 17 00:00:00 2001 From: MathieuCarriere Date: Wed, 29 Apr 2020 18:31:24 -0400 Subject: small modifs --- .../include/gudhi/Bitmap_cubical_complex_base.h | 4 +-- src/python/gudhi/cubical_complex.pyx | 42 +++++++++++----------- .../include/Persistent_cohomology_interface.h | 7 +++- 3 files changed, 29 insertions(+), 24 deletions(-) (limited to 'src') diff --git a/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_base.h b/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_base.h index 248ebdb6..eaf8a0b6 100644 --- a/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_base.h +++ b/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_base.h @@ -116,7 +116,7 @@ class Bitmap_cubical_complex_base { * be an edge or vertex of a cube. On the other hand, the output is always indicating the position of * a cube in the data structure. **/ - inline int get_top_dimensional_coface_of_a_cell(int splx); + inline int get_top_dimensional_coface_of_a_cell(size_t splx); /** * This procedure compute incidence numbers between cubes. For a cube \f$A\f$ of @@ -612,7 +612,7 @@ void Bitmap_cubical_complex_base::setup_bitmap_based_on_top_dimensional_cells } template -int Bitmap_cubical_complex_base::get_top_dimensional_coface_of_a_cell(int splx) { +int Bitmap_cubical_complex_base::get_top_dimensional_coface_of_a_cell(size_t splx) { if (this->get_dimension_of_a_cell(splx) == this->dimension()){return splx;} else{ for (auto v : this->get_coboundary_of_a_cell(splx)){ diff --git a/src/python/gudhi/cubical_complex.pyx b/src/python/gudhi/cubical_complex.pyx index 884b0664..b16a037f 100644 --- a/src/python/gudhi/cubical_complex.pyx +++ b/src/python/gudhi/cubical_complex.pyx @@ -199,28 +199,28 @@ cdef class CubicalComplex: The indices of the arrays in the list correspond to the homological dimensions, and the integers of each row in each array correspond to: (index of positive top-dimensional cell). """ + + assert self.pcohptr != NULL, "cofaces_of_persistence_pairs function requires persistence function to be launched first." + cdef vector[vector[int]] persistence_result - if self.pcohptr != NULL: - output = [[],[]] - persistence_result = self.pcohptr.cofaces_of_cubical_persistence_pairs() - pr = np.array(persistence_result) - - ess_ind = np.argwhere(pr[:,2] == -1)[:,0] - ess = pr[ess_ind] - max_h = max(ess[:,0])+1 - for h in range(max_h): - hidxs = np.argwhere(ess[:,0] == h)[:,0] - output[1].append(ess[hidxs][:,1]) - - reg_ind = np.setdiff1d(np.array(range(len(pr))), ess_ind) - reg = pr[reg_ind] - max_h = max(reg[:,0])+1 - for h in range(max_h): - hidxs = np.argwhere(reg[:,0] == h)[:,0] - output[0].append(reg[hidxs][:,1:]) - else: - print("cofaces_of_persistence_pairs function requires persistence function" - " to be launched first.") + output = [[],[]] + persistence_result = self.pcohptr.cofaces_of_cubical_persistence_pairs() + pr = np.array(persistence_result) + + ess_ind = np.argwhere(pr[:,2] == -1)[:,0] + ess = pr[ess_ind] + max_h = max(ess[:,0])+1 + for h in range(max_h): + hidxs = np.argwhere(ess[:,0] == h)[:,0] + output[1].append(ess[hidxs][:,1]) + + reg_ind = np.setdiff1d(np.array(range(len(pr))), ess_ind) + reg = pr[reg_ind] + max_h = max(reg[:,0])+1 + for h in range(max_h): + hidxs = np.argwhere(reg[:,0] == h)[:,0] + output[0].append(reg[hidxs][:,1:]) + return output def betti_numbers(self): diff --git a/src/python/include/Persistent_cohomology_interface.h b/src/python/include/Persistent_cohomology_interface.h index c4e60a27..cec18546 100644 --- a/src/python/include/Persistent_cohomology_interface.h +++ b/src/python/include/Persistent_cohomology_interface.h @@ -68,11 +68,16 @@ persistent_cohomology::Persistent_cohomology> cofaces_of_cubical_persistence_pairs() { // Warning: this function is meant to be used with CubicalComplex only!! - auto pairs = persistent_cohomology::Persistent_cohomology::get_persistent_pairs(); // Gather all top-dimensional cells and store their simplex handles -- cgit v1.2.3 From a51f4f177e29ad5b01e58c9d8dd2560fb9b4fb19 Mon Sep 17 00:00:00 2001 From: MathieuCarriere Date: Thu, 30 Apr 2020 00:52:52 -0400 Subject: int to size_t --- .../include/gudhi/Bitmap_cubical_complex_base.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) (limited to 'src') diff --git a/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_base.h b/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_base.h index eaf8a0b6..e0c567ae 100644 --- a/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_base.h +++ b/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_base.h @@ -116,7 +116,7 @@ class Bitmap_cubical_complex_base { * be an edge or vertex of a cube. On the other hand, the output is always indicating the position of * a cube in the data structure. **/ - inline int get_top_dimensional_coface_of_a_cell(size_t splx); + inline size_t get_top_dimensional_coface_of_a_cell(size_t splx); /** * This procedure compute incidence numbers between cubes. For a cube \f$A\f$ of @@ -612,7 +612,7 @@ void Bitmap_cubical_complex_base::setup_bitmap_based_on_top_dimensional_cells } template -int Bitmap_cubical_complex_base::get_top_dimensional_coface_of_a_cell(size_t splx) { +size_t Bitmap_cubical_complex_base::get_top_dimensional_coface_of_a_cell(size_t splx) { if (this->get_dimension_of_a_cell(splx) == this->dimension()){return splx;} else{ for (auto v : this->get_coboundary_of_a_cell(splx)){ -- cgit v1.2.3 From 5040c75893cb864f5e780b6644b8097f7beeb3a6 Mon Sep 17 00:00:00 2001 From: yuichi-ike Date: Mon, 11 May 2020 10:45:02 +0900 Subject: document and comments added, weights modified --- src/python/doc/rips_complex_ref.rst | 51 +++++++++++++++++++++++++++++++ src/python/gudhi/weighted_rips_complex.py | 18 ++++++----- src/python/test/test_weighted_rips.py | 2 +- 3 files changed, 63 insertions(+), 8 deletions(-) (limited to 'src') diff --git a/src/python/doc/rips_complex_ref.rst b/src/python/doc/rips_complex_ref.rst index 22b5616c..8fc7e1b0 100644 --- a/src/python/doc/rips_complex_ref.rst +++ b/src/python/doc/rips_complex_ref.rst @@ -12,3 +12,54 @@ Rips complex reference manual :show-inheritance: .. automethod:: gudhi.RipsComplex.__init__ + +====================================== +Weighted Rips complex reference manual +====================================== + +.. autoclass:: gudhi.WeightedRipsComplex + :members: + :undoc-members: + :show-inheritance: + + .. automethod:: gudhi.WeightedRipsComplex.__init__ + +Basic examples +------------- + +The following example computes the weighted Rips filtration associated with a distance matrix and weights on vertices. + +.. testcode:: + + from gudhi.weighted_rips_complex import WeightedRipsComplex + dist = [[], [1]] + weights = [1, 100] + w_rips = WeightedRipsComplex(distance_matrix=dist, weights=weights) + st = w_rips.create_simplex_tree(max_dimension=2) + print(st.get_filtration()) + +The output is: + +.. testoutput:: + + [([0], 2.0), ([1], 200.0), ([0, 1], 200.0)] + +Combining with DistanceToMeasure, one can compute the DTM-filtration of a point set, as in `this notebook `_. + +.. testcode:: + + import numpy as np + from scipy.spatial.distance import cdist + from gudhi.point_cloud.dtm import DistanceToMeasure + from gudhi.weighted_rips_complex import WeightedRipsComplex + pts = np.array([[2.0, 2.0], [0.0, 1.0], [3.0, 4.0]]) + dist = cdist(pts,pts) + dtm = DistanceToMeasure(2, q=2, metric="precomputed") + r = dtm.fit_transform(dist) + w_rips = WeightedRipsComplex(distance_matrix=dist, weights=r) + st = w_rips.create_simplex_tree(max_dimension=2) + print(st.persistence()) + +.. testoutput:: + + [(0, (3.1622776601683795, inf)), (0, (3.1622776601683795, 5.39834563766817)), (0, (3.1622776601683795, 5.39834563766817))] diff --git a/src/python/gudhi/weighted_rips_complex.py b/src/python/gudhi/weighted_rips_complex.py index 83fa82c5..7401c428 100644 --- a/src/python/gudhi/weighted_rips_complex.py +++ b/src/python/gudhi/weighted_rips_complex.py @@ -11,23 +11,26 @@ from gudhi import SimplexTree class WeightedRipsComplex: """ - Class to generate a weighted Rips complex from a distance matrix and weights on vertices. + Class to generate a weighted Rips complex from a distance matrix and weights on vertices, + in the way described in the paper 'DTM-based filtrations' https://arxiv.org/abs/1811.04757. + Remark that the filtration value of a vertex is twice of its weight for the consistency with + RipsComplex, which is different from the definition in the paper. """ def __init__(self, distance_matrix, - weights="diagonal", + weights=None, max_filtration=float('inf')): """ Args: - distance_matrix (list of list of float): distance matrix (full square or lower triangular). - weights (list of float): (one half of) weight for each vertex. + distance_matrix (Sequence[Sequence[float]]): distance matrix (full square or lower triangular). + weights (Sequence[float]): (one half of) weight for each vertex. max_filtration (float): specifies the maximal filtration value to be considered. """ self.distance_matrix = distance_matrix - if weights == "diagonal": - self.weights = [distance_matrix[i][i] for i in range(len(distance_matrix))] - else: + if weights is not None: self.weights = weights + else: + self.weights = [0] * len(distance_matrix) self.max_filtration = max_filtration def create_simplex_tree(self, max_dimension): @@ -47,6 +50,7 @@ class WeightedRipsComplex: for i in range(num_pts): for j in range(i): value = max(2*F[i], 2*F[j], dist[i][j] + F[i] + F[j]) + # max is needed when F is not 1-Lipschitz if value <= self.max_filtration: st.insert([i,j], filtration=value) diff --git a/src/python/test/test_weighted_rips.py b/src/python/test/test_weighted_rips.py index d3721115..59ec022a 100644 --- a/src/python/test/test_weighted_rips.py +++ b/src/python/test/test_weighted_rips.py @@ -51,7 +51,7 @@ def test_compatibility_with_filtered_rips(): assert st.num_vertices() == 4 def test_dtm_rips_complex(): - pts = np.array([[2.0, 2], [0, 1], [3, 4]]) + pts = np.array([[2.0, 2.0], [0.0, 1.0], [3.0, 4.0]]) dist = cdist(pts,pts) dtm = DistanceToMeasure(2, q=2, metric="precomputed") r = dtm.fit_transform(dist) -- cgit v1.2.3 From 6c17494e02721ca826750155bac14c7f91a173fa Mon Sep 17 00:00:00 2001 From: yuichi-ike Date: Tue, 12 May 2020 09:37:32 +0900 Subject: reference and comments added --- biblio/bibliography.bib | 26 ++++++++++++++++++++++++++ src/python/CMakeLists.txt | 4 +++- src/python/doc/rips_complex_ref.rst | 4 +++- src/python/gudhi/weighted_rips_complex.py | 6 +++--- src/python/test/test_weighted_rips.py | 4 ++-- 5 files changed, 37 insertions(+), 7 deletions(-) (limited to 'src') diff --git a/biblio/bibliography.bib b/biblio/bibliography.bib index 99a15c5e..f405b9bb 100644 --- a/biblio/bibliography.bib +++ b/biblio/bibliography.bib @@ -1247,3 +1247,29 @@ year = "2011" year={2014}, publisher={Springer} } + +@inproceedings{dtmfiltrations, + author = {Hirokazu Anai and + Fr{\'{e}}d{\'{e}}ric Chazal and + Marc Glisse and + Yuichi Ike and + Hiroya Inakoshi and + Rapha{\"{e}}l Tinarrage and + Yuhei Umeda}, + editor = {Gill Barequet and + Yusu Wang}, + title = {DTM-Based Filtrations}, + booktitle = {35th International Symposium on Computational Geometry, SoCG 2019, + June 18-21, 2019, Portland, Oregon, {USA}}, + series = {LIPIcs}, + volume = {129}, + pages = {58:1--58:15}, + publisher = {Schloss Dagstuhl - Leibniz-Zentrum f{\"{u}}r Informatik}, + year = {2019}, + url = {https://doi.org/10.4230/LIPIcs.SoCG.2019.58}, + doi = {10.4230/LIPIcs.SoCG.2019.58}, + timestamp = {Tue, 11 Feb 2020 15:52:14 +0100}, + biburl = {https://dblp.org/rec/conf/compgeom/AnaiCGIITU19.bib}, + bibsource = {dblp computer science bibliography, https://dblp.org} +} + diff --git a/src/python/CMakeLists.txt b/src/python/CMakeLists.txt index adf4923b..0aa55467 100644 --- a/src/python/CMakeLists.txt +++ b/src/python/CMakeLists.txt @@ -487,7 +487,9 @@ if(PYTHONINTERP_FOUND) endif() # Weighted Rips - add_gudhi_py_test(test_weighted_rips) + if(SCIPY_FOUND) + add_gudhi_py_test(test_weighted_rips) + endif() # Set missing or not modules set(GUDHI_MODULES ${GUDHI_MODULES} "python" CACHE INTERNAL "GUDHI_MODULES") diff --git a/src/python/doc/rips_complex_ref.rst b/src/python/doc/rips_complex_ref.rst index 8fc7e1b0..3c25564a 100644 --- a/src/python/doc/rips_complex_ref.rst +++ b/src/python/doc/rips_complex_ref.rst @@ -25,7 +25,7 @@ Weighted Rips complex reference manual .. automethod:: gudhi.WeightedRipsComplex.__init__ Basic examples -------------- +-------------- The following example computes the weighted Rips filtration associated with a distance matrix and weights on vertices. @@ -60,6 +60,8 @@ Combining with DistanceToMeasure, one can compute the DTM-filtration of a point st = w_rips.create_simplex_tree(max_dimension=2) print(st.persistence()) +The output is: + .. testoutput:: [(0, (3.1622776601683795, inf)), (0, (3.1622776601683795, 5.39834563766817)), (0, (3.1622776601683795, 5.39834563766817))] diff --git a/src/python/gudhi/weighted_rips_complex.py b/src/python/gudhi/weighted_rips_complex.py index 7401c428..bccac1ff 100644 --- a/src/python/gudhi/weighted_rips_complex.py +++ b/src/python/gudhi/weighted_rips_complex.py @@ -12,9 +12,9 @@ from gudhi import SimplexTree class WeightedRipsComplex: """ Class to generate a weighted Rips complex from a distance matrix and weights on vertices, - in the way described in the paper 'DTM-based filtrations' https://arxiv.org/abs/1811.04757. - Remark that the filtration value of a vertex is twice of its weight for the consistency with - RipsComplex, which is different from the definition in the paper. + in the way described in :cite:`dtmfiltrations`. + Remark that all the filtration values of vertices are twice of the given weights for the consistency + with RipsComplex, which is different from the definition in the paper. """ def __init__(self, distance_matrix, diff --git a/src/python/test/test_weighted_rips.py b/src/python/test/test_weighted_rips.py index 59ec022a..7ef48333 100644 --- a/src/python/test/test_weighted_rips.py +++ b/src/python/test/test_weighted_rips.py @@ -35,8 +35,8 @@ def test_compatibility_with_rips(): ([0, 2], 1.0), ([1, 3], 1.0), ([2, 3], 1.0), - ([1, 2], 1.4142135623730951), - ([0, 3], 1.4142135623730951), + ([1, 2], sqrt(2)), + ([0, 3], sqrt(2)), ] def test_compatibility_with_filtered_rips(): -- cgit v1.2.3 From a9c1e13e7f994e5c8d9f1c3d0311a5815df1e67d Mon Sep 17 00:00:00 2001 From: yuichi-ike Date: Tue, 12 May 2020 11:10:16 +0900 Subject: document fixed --- src/python/doc/rips_complex_ref.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) (limited to 'src') diff --git a/src/python/doc/rips_complex_ref.rst b/src/python/doc/rips_complex_ref.rst index 3c25564a..8946d156 100644 --- a/src/python/doc/rips_complex_ref.rst +++ b/src/python/doc/rips_complex_ref.rst @@ -22,7 +22,7 @@ Weighted Rips complex reference manual :undoc-members: :show-inheritance: - .. automethod:: gudhi.WeightedRipsComplex.__init__ + .. automethod:: gudhi.weighted_rips_complex.WeightedRipsComplex.__init__ Basic examples -------------- -- cgit v1.2.3 From 23547c0cbbe9e42b4dfadec3a116751302fd19ab Mon Sep 17 00:00:00 2001 From: yuichi-ike Date: Tue, 12 May 2020 11:41:03 +0900 Subject: document fixed --- src/python/doc/rips_complex_ref.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) (limited to 'src') diff --git a/src/python/doc/rips_complex_ref.rst b/src/python/doc/rips_complex_ref.rst index 8946d156..1f73f95b 100644 --- a/src/python/doc/rips_complex_ref.rst +++ b/src/python/doc/rips_complex_ref.rst @@ -17,7 +17,7 @@ Rips complex reference manual Weighted Rips complex reference manual ====================================== -.. autoclass:: gudhi.WeightedRipsComplex +.. autoclass:: gudhi.weighted_rips_complex.WeightedRipsComplex :members: :undoc-members: :show-inheritance: -- cgit v1.2.3 From 2c4049895bb2844c2ad1b43b9df51ad5b259fc39 Mon Sep 17 00:00:00 2001 From: yuichi-ike Date: Tue, 12 May 2020 13:09:40 +0900 Subject: a test in a document fixed --- src/python/doc/rips_complex_ref.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) (limited to 'src') diff --git a/src/python/doc/rips_complex_ref.rst b/src/python/doc/rips_complex_ref.rst index 1f73f95b..a5b4ffed 100644 --- a/src/python/doc/rips_complex_ref.rst +++ b/src/python/doc/rips_complex_ref.rst @@ -36,7 +36,7 @@ The following example computes the weighted Rips filtration associated with a di weights = [1, 100] w_rips = WeightedRipsComplex(distance_matrix=dist, weights=weights) st = w_rips.create_simplex_tree(max_dimension=2) - print(st.get_filtration()) + print(list(st.get_filtration())) The output is: -- cgit v1.2.3 From c60caee5623d0b1ef55e7b2a5854604080419df1 Mon Sep 17 00:00:00 2001 From: yuichi-ike Date: Tue, 12 May 2020 15:06:55 +0900 Subject: comment modified --- src/python/gudhi/weighted_rips_complex.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) (limited to 'src') diff --git a/src/python/gudhi/weighted_rips_complex.py b/src/python/gudhi/weighted_rips_complex.py index bccac1ff..0541572b 100644 --- a/src/python/gudhi/weighted_rips_complex.py +++ b/src/python/gudhi/weighted_rips_complex.py @@ -13,8 +13,8 @@ class WeightedRipsComplex: """ Class to generate a weighted Rips complex from a distance matrix and weights on vertices, in the way described in :cite:`dtmfiltrations`. - Remark that all the filtration values of vertices are twice of the given weights for the consistency - with RipsComplex, which is different from the definition in the paper. + Remark that all the filtration values are doubled compared to the definition in the paper + for the consistency with RipsComplex. """ def __init__(self, distance_matrix, -- cgit v1.2.3 From fd7112b7e665d495543d9647f675a14f75061bbf Mon Sep 17 00:00:00 2001 From: yuichi-ike Date: Wed, 13 May 2020 09:54:47 +0900 Subject: documents modified --- src/python/doc/rips_complex_ref.rst | 42 ------------------------------- src/python/doc/rips_complex_sum.inc | 3 +++ src/python/doc/rips_complex_user.rst | 48 ++++++++++++++++++++++++++++++++++++ 3 files changed, 51 insertions(+), 42 deletions(-) (limited to 'src') diff --git a/src/python/doc/rips_complex_ref.rst b/src/python/doc/rips_complex_ref.rst index a5b4ffed..9ae3c49c 100644 --- a/src/python/doc/rips_complex_ref.rst +++ b/src/python/doc/rips_complex_ref.rst @@ -23,45 +23,3 @@ Weighted Rips complex reference manual :show-inheritance: .. automethod:: gudhi.weighted_rips_complex.WeightedRipsComplex.__init__ - -Basic examples --------------- - -The following example computes the weighted Rips filtration associated with a distance matrix and weights on vertices. - -.. testcode:: - - from gudhi.weighted_rips_complex import WeightedRipsComplex - dist = [[], [1]] - weights = [1, 100] - w_rips = WeightedRipsComplex(distance_matrix=dist, weights=weights) - st = w_rips.create_simplex_tree(max_dimension=2) - print(list(st.get_filtration())) - -The output is: - -.. testoutput:: - - [([0], 2.0), ([1], 200.0), ([0, 1], 200.0)] - -Combining with DistanceToMeasure, one can compute the DTM-filtration of a point set, as in `this notebook `_. - -.. testcode:: - - import numpy as np - from scipy.spatial.distance import cdist - from gudhi.point_cloud.dtm import DistanceToMeasure - from gudhi.weighted_rips_complex import WeightedRipsComplex - pts = np.array([[2.0, 2.0], [0.0, 1.0], [3.0, 4.0]]) - dist = cdist(pts,pts) - dtm = DistanceToMeasure(2, q=2, metric="precomputed") - r = dtm.fit_transform(dist) - w_rips = WeightedRipsComplex(distance_matrix=dist, weights=r) - st = w_rips.create_simplex_tree(max_dimension=2) - print(st.persistence()) - -The output is: - -.. testoutput:: - - [(0, (3.1622776601683795, inf)), (0, (3.1622776601683795, 5.39834563766817)), (0, (3.1622776601683795, 5.39834563766817))] diff --git a/src/python/doc/rips_complex_sum.inc b/src/python/doc/rips_complex_sum.inc index 6feb74cd..f7580714 100644 --- a/src/python/doc/rips_complex_sum.inc +++ b/src/python/doc/rips_complex_sum.inc @@ -11,6 +11,9 @@ | | | | | | This complex can be built from a point cloud and a distance function, | | | | or from a distance matrix. | | + | | | | + | | Weighted Rips complex constructs a simplicial complex from a distance | | + | | matrix and weights on vertices. | | +----------------------------------------------------------------+------------------------------------------------------------------------+----------------------------------------------------------------------+ | * :doc:`rips_complex_user` | * :doc:`rips_complex_ref` | +----------------------------------------------------------------+-----------------------------------------------------------------------------------------------------------------------------------------------+ diff --git a/src/python/doc/rips_complex_user.rst b/src/python/doc/rips_complex_user.rst index 8efb12e6..adb002a8 100644 --- a/src/python/doc/rips_complex_user.rst +++ b/src/python/doc/rips_complex_user.rst @@ -347,3 +347,51 @@ until dimension 1 - one skeleton graph in other words), the output is: points in the persistence diagram will be under the diagonal, and bottleneck distance and persistence graphical tool will not work properly, this is a known issue. + +Weighted Rips Complex +--------------------- + +Example from a distance matrix and weights +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +The following example computes the weighted Rips filtration associated with a distance matrix and weights on vertices. + +.. testcode:: + + from gudhi.weighted_rips_complex import WeightedRipsComplex + dist = [[], [1]] + weights = [1, 100] + w_rips = WeightedRipsComplex(distance_matrix=dist, weights=weights) + st = w_rips.create_simplex_tree(max_dimension=2) + print(list(st.get_filtration())) + +The output is: + +.. testoutput:: + + [([0], 2.0), ([1], 200.0), ([0, 1], 200.0)] + +Example from a point cloud combined with DistanceToMeasure +---------------------------------------------------------- + +Combining with DistanceToMeasure, one can compute the DTM-filtration of a point set, as in `this notebook `_. + +.. testcode:: + + import numpy as np + from scipy.spatial.distance import cdist + from gudhi.point_cloud.dtm import DistanceToMeasure + from gudhi.weighted_rips_complex import WeightedRipsComplex + pts = np.array([[2.0, 2.0], [0.0, 1.0], [3.0, 4.0]]) + dist = cdist(pts,pts) + dtm = DistanceToMeasure(2, q=2, metric="precomputed") + r = dtm.fit_transform(dist) + w_rips = WeightedRipsComplex(distance_matrix=dist, weights=r) + st = w_rips.create_simplex_tree(max_dimension=2) + print(st.persistence()) + +The output is: + +.. testoutput:: + + [(0, (3.1622776601683795, inf)), (0, (3.1622776601683795, 5.39834563766817)), (0, (3.1622776601683795, 5.39834563766817))] -- cgit v1.2.3 From 7b4ffb762edae9036cbec12b34eeb64f2cffd0e7 Mon Sep 17 00:00:00 2001 From: Marc Glisse Date: Wed, 13 May 2020 18:12:58 +0200 Subject: Rephrase comment about cubes --- .../include/gudhi/Bitmap_cubical_complex_base.h | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) (limited to 'src') diff --git a/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_base.h b/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_base.h index e0c567ae..99487dc3 100644 --- a/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_base.h +++ b/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_base.h @@ -112,9 +112,9 @@ class Bitmap_cubical_complex_base { /** * This function finds a top-dimensional cell that is incident to the input cell and has * the same filtration value. In case several cells are suitable, an arbitrary one is - * returned. Note that the input parameter is not necessarily a cube, it might also - * be an edge or vertex of a cube. On the other hand, the output is always indicating the position of - * a cube in the data structure. + * returned. Note that the input parameter can be a cell of any dimension (vertex, edge, etc). + * On the other hand, the output is always indicating the position of + * a top-dimensional cube in the data structure. **/ inline size_t get_top_dimensional_coface_of_a_cell(size_t splx); -- cgit v1.2.3 From 5c3e042628b7db2b82d92f644f7ab0fc409a357b Mon Sep 17 00:00:00 2001 From: Marc Glisse Date: Wed, 13 May 2020 18:45:28 +0200 Subject: BOOST_UNREACHABLE_RETURN + comment --- .../include/gudhi/Bitmap_cubical_complex_base.h | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) (limited to 'src') diff --git a/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_base.h b/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_base.h index 99487dc3..5927bbec 100644 --- a/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_base.h +++ b/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_base.h @@ -13,6 +13,8 @@ #include +#include + #include #include #include @@ -115,6 +117,7 @@ class Bitmap_cubical_complex_base { * returned. Note that the input parameter can be a cell of any dimension (vertex, edge, etc). * On the other hand, the output is always indicating the position of * a top-dimensional cube in the data structure. + * \pre The filtration values are assigned as per `impose_lower_star_filtration()`. **/ inline size_t get_top_dimensional_coface_of_a_cell(size_t splx); @@ -621,7 +624,7 @@ size_t Bitmap_cubical_complex_base::get_top_dimensional_coface_of_a_cell(size } } } - return splx; + BOOST_UNREACHABLE_RETURN(-2); } template -- cgit v1.2.3 From b2118cde83056b43cea095f5208d37744c9f088f Mon Sep 17 00:00:00 2001 From: Marc Glisse Date: Wed, 13 May 2020 18:51:16 +0200 Subject: compute_persistence --- src/python/gudhi/cubical_complex.pyx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) (limited to 'src') diff --git a/src/python/gudhi/cubical_complex.pyx b/src/python/gudhi/cubical_complex.pyx index b16a037f..9ebd0b30 100644 --- a/src/python/gudhi/cubical_complex.pyx +++ b/src/python/gudhi/cubical_complex.pyx @@ -200,7 +200,7 @@ cdef class CubicalComplex: integers of each row in each array correspond to: (index of positive top-dimensional cell). """ - assert self.pcohptr != NULL, "cofaces_of_persistence_pairs function requires persistence function to be launched first." + assert self.pcohptr != NULL, "compute_persistence() must be called before cofaces_of_persistence_pairs()" cdef vector[vector[int]] persistence_result output = [[],[]] -- cgit v1.2.3 From 7bbc1ae35d492123c517a54a9595188938e52dff Mon Sep 17 00:00:00 2001 From: Marc Glisse Date: Wed, 13 May 2020 19:32:21 +0200 Subject: More size_t --- .../include/Persistent_cohomology_interface.h | 28 +++++++++++----------- 1 file changed, 14 insertions(+), 14 deletions(-) (limited to 'src') diff --git a/src/python/include/Persistent_cohomology_interface.h b/src/python/include/Persistent_cohomology_interface.h index cec18546..e5a3dfba 100644 --- a/src/python/include/Persistent_cohomology_interface.h +++ b/src/python/include/Persistent_cohomology_interface.h @@ -13,6 +13,7 @@ #include +#include #include #include // for std::pair #include // for sort @@ -81,32 +82,31 @@ persistent_cohomology::Persistent_cohomology::get_persistent_pairs(); // Gather all top-dimensional cells and store their simplex handles - std::vector max_splx; for (auto splx : stptr_->top_dimensional_cells_range()){ - max_splx.push_back(splx); - } + std::vector max_splx; + for (auto splx : stptr_->top_dimensional_cells_range()) + max_splx.push_back(splx); // Sort these simplex handles and compute the ordering function - // This function allows to go directly from the simplex handle to the position of the corresponding top-dimensional cell in the input data - std::unordered_map order; - //std::sort(max_splx.begin(), max_splx.end()); + // This function allows to go directly from the simplex handle to the position of the corresponding top-dimensional cell in the input data + std::unordered_map order; + //std::sort(max_splx.begin(), max_splx.end()); for (unsigned int i = 0; i < max_splx.size(); i++) order.emplace(max_splx[i], i); std::vector> persistence_pairs; for (auto pair : pairs) { int h = stptr_->dimension(get<0>(pair)); // Recursively get the top-dimensional cell / coface associated to the persistence generator - int face0 = stptr_->get_top_dimensional_coface_of_a_cell(get<0>(pair)); + std::size_t face0 = stptr_->get_top_dimensional_coface_of_a_cell(get<0>(pair)); // Retrieve the index of the corresponding top-dimensional cell in the input data int splx0 = order[face0]; int splx1 = -1; - if (isfinite(stptr_->filtration(get<1>(pair)))){ - // Recursively get the top-dimensional cell / coface associated to the persistence generator - int face1 = stptr_->get_top_dimensional_coface_of_a_cell(get<1>(pair)); - // Retrieve the index of the corresponding top-dimensional cell in the input data - splx1 = order[face1]; + if (get<1>(pair) != stptr_->null_simplex()){ + // Recursively get the top-dimensional cell / coface associated to the persistence generator + std::size_t face1 = stptr_->get_top_dimensional_coface_of_a_cell(get<1>(pair)); + // Retrieve the index of the corresponding top-dimensional cell in the input data + splx1 = order[face1]; } - std::vector vect{ h, splx0, splx1}; - persistence_pairs.push_back(vect); + persistence_pairs.push_back({ h, splx0, splx1 }); } return persistence_pairs; } -- cgit v1.2.3 From b0ae08e93fdba8a1faec56c2230b6f542653c49e Mon Sep 17 00:00:00 2001 From: Marc Glisse Date: Wed, 13 May 2020 20:17:26 +0200 Subject: Trailing whitespace --- .../include/gudhi/Bitmap_cubical_complex_base.h | 8 ++--- src/python/gudhi/cubical_complex.pyx | 34 +++++++++++----------- src/python/gudhi/periodic_cubical_complex.pyx | 34 +++++++++++----------- 3 files changed, 38 insertions(+), 38 deletions(-) (limited to 'src') diff --git a/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_base.h b/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_base.h index 5927bbec..58d9208d 100644 --- a/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_base.h +++ b/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_base.h @@ -112,8 +112,8 @@ class Bitmap_cubical_complex_base { virtual inline std::vector get_coboundary_of_a_cell(std::size_t cell) const; /** - * This function finds a top-dimensional cell that is incident to the input cell and has - * the same filtration value. In case several cells are suitable, an arbitrary one is + * This function finds a top-dimensional cell that is incident to the input cell and has + * the same filtration value. In case several cells are suitable, an arbitrary one is * returned. Note that the input parameter can be a cell of any dimension (vertex, edge, etc). * On the other hand, the output is always indicating the position of * a top-dimensional cube in the data structure. @@ -617,12 +617,12 @@ void Bitmap_cubical_complex_base::setup_bitmap_based_on_top_dimensional_cells template size_t Bitmap_cubical_complex_base::get_top_dimensional_coface_of_a_cell(size_t splx) { if (this->get_dimension_of_a_cell(splx) == this->dimension()){return splx;} - else{ + else{ for (auto v : this->get_coboundary_of_a_cell(splx)){ if(this->get_cell_data(v) == this->get_cell_data(splx)){ return this->get_top_dimensional_coface_of_a_cell(v); } - } + } } BOOST_UNREACHABLE_RETURN(-2); } diff --git a/src/python/gudhi/cubical_complex.pyx b/src/python/gudhi/cubical_complex.pyx index 9ebd0b30..ca979eda 100644 --- a/src/python/gudhi/cubical_complex.pyx +++ b/src/python/gudhi/cubical_complex.pyx @@ -172,31 +172,31 @@ cdef class CubicalComplex: return self.pcohptr.get_persistence() def cofaces_of_persistence_pairs(self): - """A persistence interval is described by a pair of cells, one that creates the - feature and one that kills it. The filtration values of those 2 cells give coordinates - for a point in a persistence diagram, or a bar in a barcode. Structurally, in the - cubical complexes provided here, the filtration value of any cell is the minimum of the - filtration values of the maximal cells that contain it. Connecting persistence diagram - coordinates to the corresponding value in the input (i.e. the filtration values of + """A persistence interval is described by a pair of cells, one that creates the + feature and one that kills it. The filtration values of those 2 cells give coordinates + for a point in a persistence diagram, or a bar in a barcode. Structurally, in the + cubical complexes provided here, the filtration value of any cell is the minimum of the + filtration values of the maximal cells that contain it. Connecting persistence diagram + coordinates to the corresponding value in the input (i.e. the filtration values of the top-dimensional cells) is useful for differentiation purposes. - This function returns a list of pairs of top-dimensional cells corresponding to - the persistence birth and death cells of the filtration. The cells are represented by - their indices in the input list of top-dimensional cells (and not their indices in the - internal datastructure that includes non-maximal cells). Note that when two adjacent + This function returns a list of pairs of top-dimensional cells corresponding to + the persistence birth and death cells of the filtration. The cells are represented by + their indices in the input list of top-dimensional cells (and not their indices in the + internal datastructure that includes non-maximal cells). Note that when two adjacent top-dimensional cells have the same filtration value, we arbitrarily return one of the two when calling the function on one of their common faces. - :returns: The top-dimensional cells/cofaces of the positive and negative cells, + :returns: The top-dimensional cells/cofaces of the positive and negative cells, together with the corresponding homological dimension, in two lists of numpy arrays of integers. - The first list contains the regular persistence pairs, grouped by dimension. + The first list contains the regular persistence pairs, grouped by dimension. It contains numpy arrays of shape [number_of_persistence_points, 2]. - The indices of the arrays in the list correspond to the homological dimensions, and the - integers of each row in each array correspond to: (index of positive top-dimensional cell, - index of negative top-dimensional cell). - The second list contains the essential features, grouped by dimension. + The indices of the arrays in the list correspond to the homological dimensions, and the + integers of each row in each array correspond to: (index of positive top-dimensional cell, + index of negative top-dimensional cell). + The second list contains the essential features, grouped by dimension. It contains numpy arrays of shape [number_of_persistence_points, 1]. - The indices of the arrays in the list correspond to the homological dimensions, and the + The indices of the arrays in the list correspond to the homological dimensions, and the integers of each row in each array correspond to: (index of positive top-dimensional cell). """ diff --git a/src/python/gudhi/periodic_cubical_complex.pyx b/src/python/gudhi/periodic_cubical_complex.pyx index 3cf2ff01..06309772 100644 --- a/src/python/gudhi/periodic_cubical_complex.pyx +++ b/src/python/gudhi/periodic_cubical_complex.pyx @@ -177,31 +177,31 @@ cdef class PeriodicCubicalComplex: return self.pcohptr.get_persistence() def cofaces_of_persistence_pairs(self): - """A persistence interval is described by a pair of cells, one that creates the - feature and one that kills it. The filtration values of those 2 cells give coordinates - for a point in a persistence diagram, or a bar in a barcode. Structurally, in the - cubical complexes provided here, the filtration value of any cell is the minimum of the - filtration values of the maximal cells that contain it. Connecting persistence diagram - coordinates to the corresponding value in the input (i.e. the filtration values of + """A persistence interval is described by a pair of cells, one that creates the + feature and one that kills it. The filtration values of those 2 cells give coordinates + for a point in a persistence diagram, or a bar in a barcode. Structurally, in the + cubical complexes provided here, the filtration value of any cell is the minimum of the + filtration values of the maximal cells that contain it. Connecting persistence diagram + coordinates to the corresponding value in the input (i.e. the filtration values of the top-dimensional cells) is useful for differentiation purposes. - This function returns a list of pairs of top-dimensional cells corresponding to - the persistence birth and death cells of the filtration. The cells are represented by - their indices in the input list of top-dimensional cells (and not their indices in the - internal datastructure that includes non-maximal cells). Note that when two adjacent + This function returns a list of pairs of top-dimensional cells corresponding to + the persistence birth and death cells of the filtration. The cells are represented by + their indices in the input list of top-dimensional cells (and not their indices in the + internal datastructure that includes non-maximal cells). Note that when two adjacent top-dimensional cells have the same filtration value, we arbitrarily return one of the two when calling the function on one of their common faces. - :returns: The top-dimensional cells/cofaces of the positive and negative cells, + :returns: The top-dimensional cells/cofaces of the positive and negative cells, together with the corresponding homological dimension, in two lists of numpy arrays of integers. - The first list contains the regular persistence pairs, grouped by dimension. + The first list contains the regular persistence pairs, grouped by dimension. It contains numpy arrays of shape [number_of_persistence_points, 2]. - The indices of the arrays in the list correspond to the homological dimensions, and the - integers of each row in each array correspond to: (index of positive top-dimensional cell, - index of negative top-dimensional cell). - The second list contains the essential features, grouped by dimension. + The indices of the arrays in the list correspond to the homological dimensions, and the + integers of each row in each array correspond to: (index of positive top-dimensional cell, + index of negative top-dimensional cell). + The second list contains the essential features, grouped by dimension. It contains numpy arrays of shape [number_of_persistence_points, 1]. - The indices of the arrays in the list correspond to the homological dimensions, and the + The indices of the arrays in the list correspond to the homological dimensions, and the integers of each row in each array correspond to: (index of positive top-dimensional cell). """ cdef vector[vector[int]] persistence_result -- cgit v1.2.3 From 4d27d32308f94e63d76bbd5564b8837b94b24339 Mon Sep 17 00:00:00 2001 From: yuichi-ike Date: Thu, 14 May 2020 17:56:10 +0900 Subject: document modified --- src/python/doc/rips_complex_ref.rst | 2 ++ src/python/doc/rips_complex_user.rst | 5 ++++- 2 files changed, 6 insertions(+), 1 deletion(-) (limited to 'src') diff --git a/src/python/doc/rips_complex_ref.rst b/src/python/doc/rips_complex_ref.rst index 9ae3c49c..5f3e46c1 100644 --- a/src/python/doc/rips_complex_ref.rst +++ b/src/python/doc/rips_complex_ref.rst @@ -13,6 +13,8 @@ Rips complex reference manual .. automethod:: gudhi.RipsComplex.__init__ +.. _weighted-rips-complex-reference-manual: + ====================================== Weighted Rips complex reference manual ====================================== diff --git a/src/python/doc/rips_complex_user.rst b/src/python/doc/rips_complex_user.rst index adb002a8..819568be 100644 --- a/src/python/doc/rips_complex_user.rst +++ b/src/python/doc/rips_complex_user.rst @@ -351,6 +351,9 @@ until dimension 1 - one skeleton graph in other words), the output is: Weighted Rips Complex --------------------- +`WeightedRipsComplex `_ builds a simplicial complex from a distance matrix and weights on vertices. + + Example from a distance matrix and weights ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ @@ -372,7 +375,7 @@ The output is: [([0], 2.0), ([1], 200.0), ([0, 1], 200.0)] Example from a point cloud combined with DistanceToMeasure ----------------------------------------------------------- +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Combining with DistanceToMeasure, one can compute the DTM-filtration of a point set, as in `this notebook `_. -- cgit v1.2.3 From 84b823b6436746a06cb8323fecd7b1f38d7ba244 Mon Sep 17 00:00:00 2001 From: Marc Glisse Date: Sat, 16 May 2020 13:30:19 +0200 Subject: Minimal nogil support for cubical complexes --- src/python/gudhi/cubical_complex.pyx | 29 +++++++++++++++------------ src/python/gudhi/periodic_cubical_complex.pyx | 29 +++++++++++++++------------ 2 files changed, 32 insertions(+), 26 deletions(-) (limited to 'src') diff --git a/src/python/gudhi/cubical_complex.pyx b/src/python/gudhi/cubical_complex.pyx index ca979eda..308b5099 100644 --- a/src/python/gudhi/cubical_complex.pyx +++ b/src/python/gudhi/cubical_complex.pyx @@ -27,20 +27,20 @@ __license__ = "MIT" cdef extern from "Cubical_complex_interface.h" namespace "Gudhi": cdef cppclass Bitmap_cubical_complex_base_interface "Gudhi::Cubical_complex::Cubical_complex_interface<>": - Bitmap_cubical_complex_base_interface(vector[unsigned] dimensions, vector[double] top_dimensional_cells) - Bitmap_cubical_complex_base_interface(string perseus_file) - int num_simplices() - int dimension() + Bitmap_cubical_complex_base_interface(vector[unsigned] dimensions, vector[double] top_dimensional_cells) nogil + Bitmap_cubical_complex_base_interface(string perseus_file) nogil + int num_simplices() nogil + int dimension() nogil cdef extern from "Persistent_cohomology_interface.h" namespace "Gudhi": cdef cppclass Cubical_complex_persistence_interface "Gudhi::Persistent_cohomology_interface>": - Cubical_complex_persistence_interface(Bitmap_cubical_complex_base_interface * st, bool persistence_dim_max) - void compute_persistence(int homology_coeff_field, double min_persistence) - vector[pair[int, pair[double, double]]] get_persistence() - vector[vector[int]] cofaces_of_cubical_persistence_pairs() - vector[int] betti_numbers() - vector[int] persistent_betti_numbers(double from_value, double to_value) - vector[pair[double,double]] intervals_in_dimension(int dimension) + Cubical_complex_persistence_interface(Bitmap_cubical_complex_base_interface * st, bool persistence_dim_max) nogil + void compute_persistence(int homology_coeff_field, double min_persistence) nogil + vector[pair[int, pair[double, double]]] get_persistence() nogil + vector[vector[int]] cofaces_of_cubical_persistence_pairs() nogil + vector[int] betti_numbers() nogil + vector[int] persistent_betti_numbers(double from_value, double to_value) nogil + vector[pair[double,double]] intervals_in_dimension(int dimension) nogil # CubicalComplex python interface cdef class CubicalComplex: @@ -151,8 +151,11 @@ cdef class CubicalComplex: if self.pcohptr != NULL: del self.pcohptr assert self.__is_defined() - self.pcohptr = new Cubical_complex_persistence_interface(self.thisptr, True) - self.pcohptr.compute_persistence(homology_coeff_field, min_persistence) + cdef int field = homology_coeff_field + cdef double minp = min_persistence + with nogil: + self.pcohptr = new Cubical_complex_persistence_interface(self.thisptr, 1) + self.pcohptr.compute_persistence(field, minp) def persistence(self, homology_coeff_field=11, min_persistence=0): """This function computes and returns the persistence of the complex. diff --git a/src/python/gudhi/periodic_cubical_complex.pyx b/src/python/gudhi/periodic_cubical_complex.pyx index 06309772..dcca7b63 100644 --- a/src/python/gudhi/periodic_cubical_complex.pyx +++ b/src/python/gudhi/periodic_cubical_complex.pyx @@ -24,20 +24,20 @@ __license__ = "MIT" cdef extern from "Cubical_complex_interface.h" namespace "Gudhi": cdef cppclass Periodic_cubical_complex_base_interface "Gudhi::Cubical_complex::Cubical_complex_interface>": - Periodic_cubical_complex_base_interface(vector[unsigned] dimensions, vector[double] top_dimensional_cells, vector[bool] periodic_dimensions) - Periodic_cubical_complex_base_interface(string perseus_file) - int num_simplices() - int dimension() + Periodic_cubical_complex_base_interface(vector[unsigned] dimensions, vector[double] top_dimensional_cells, vector[bool] periodic_dimensions) nogil + Periodic_cubical_complex_base_interface(string perseus_file) nogil + int num_simplices() nogil + int dimension() nogil cdef extern from "Persistent_cohomology_interface.h" namespace "Gudhi": cdef cppclass Periodic_cubical_complex_persistence_interface "Gudhi::Persistent_cohomology_interface>>": - Periodic_cubical_complex_persistence_interface(Periodic_cubical_complex_base_interface * st, bool persistence_dim_max) - void compute_persistence(int homology_coeff_field, double min_persistence) - vector[pair[int, pair[double, double]]] get_persistence() - vector[vector[int]] cofaces_of_cubical_persistence_pairs() - vector[int] betti_numbers() - vector[int] persistent_betti_numbers(double from_value, double to_value) - vector[pair[double,double]] intervals_in_dimension(int dimension) + Periodic_cubical_complex_persistence_interface(Periodic_cubical_complex_base_interface * st, bool persistence_dim_max) nogil + void compute_persistence(int homology_coeff_field, double min_persistence) nogil + vector[pair[int, pair[double, double]]] get_persistence() nogil + vector[vector[int]] cofaces_of_cubical_persistence_pairs() nogil + vector[int] betti_numbers() nogil + vector[int] persistent_betti_numbers(double from_value, double to_value) nogil + vector[pair[double,double]] intervals_in_dimension(int dimension) nogil # PeriodicCubicalComplex python interface cdef class PeriodicCubicalComplex: @@ -156,8 +156,11 @@ cdef class PeriodicCubicalComplex: if self.pcohptr != NULL: del self.pcohptr assert self.__is_defined() - self.pcohptr = new Periodic_cubical_complex_persistence_interface(self.thisptr, True) - self.pcohptr.compute_persistence(homology_coeff_field, min_persistence) + cdef int field = homology_coeff_field + cdef double minp = min_persistence + with nogil: + self.pcohptr = new Periodic_cubical_complex_persistence_interface(self.thisptr, 1) + self.pcohptr.compute_persistence(field, minp) def persistence(self, homology_coeff_field=11, min_persistence=0): """This function computes and returns the persistence of the complex. -- cgit v1.2.3 From 207050fb1f5af375a98c70dbd5fc22149d6f6e22 Mon Sep 17 00:00:00 2001 From: Marc Glisse Date: Sat, 16 May 2020 14:08:23 +0200 Subject: nogil for cubical constructor There may be some extra copying until cython3, but it is probably not that bad. --- src/python/gudhi/cubical_complex.pyx | 14 +++++++++++--- src/python/gudhi/periodic_cubical_complex.pyx | 18 +++++++++++------- 2 files changed, 22 insertions(+), 10 deletions(-) (limited to 'src') diff --git a/src/python/gudhi/cubical_complex.pyx b/src/python/gudhi/cubical_complex.pyx index 308b5099..0068f2ff 100644 --- a/src/python/gudhi/cubical_complex.pyx +++ b/src/python/gudhi/cubical_complex.pyx @@ -80,7 +80,7 @@ cdef class CubicalComplex: perseus_file=''): if ((dimensions is not None) and (top_dimensional_cells is not None) and (perseus_file == '')): - self.thisptr = new Bitmap_cubical_complex_base_interface(dimensions, top_dimensional_cells) + self._construct_from_cells(dimensions, top_dimensional_cells) elif ((dimensions is None) and (top_dimensional_cells is not None) and (perseus_file == '')): top_dimensional_cells = np.array(top_dimensional_cells, @@ -88,11 +88,11 @@ cdef class CubicalComplex: order = 'F') dimensions = top_dimensional_cells.shape top_dimensional_cells = top_dimensional_cells.ravel(order='F') - self.thisptr = new Bitmap_cubical_complex_base_interface(dimensions, top_dimensional_cells) + self._construct_from_cells(dimensions, top_dimensional_cells) elif ((dimensions is None) and (top_dimensional_cells is None) and (perseus_file != '')): if os.path.isfile(perseus_file): - self.thisptr = new Bitmap_cubical_complex_base_interface(perseus_file.encode('utf-8')) + self._construct_from_file(perseus_file.encode('utf-8')) else: raise FileNotFoundError(errno.ENOENT, os.strerror(errno.ENOENT), perseus_file) @@ -101,6 +101,14 @@ cdef class CubicalComplex: "top_dimensional_cells or from a Perseus-style file name.", file=sys.stderr) + def _construct_from_cells(self, vector[unsigned] dimensions, vector[double] top_dimensional_cells): + with nogil: + self.thisptr = new Bitmap_cubical_complex_base_interface(dimensions, top_dimensional_cells) + + def _construct_from_file(self, string filename): + with nogil: + self.thisptr = new Bitmap_cubical_complex_base_interface(filename) + def __dealloc__(self): if self.thisptr != NULL: del self.thisptr diff --git a/src/python/gudhi/periodic_cubical_complex.pyx b/src/python/gudhi/periodic_cubical_complex.pyx index dcca7b63..11e1766c 100644 --- a/src/python/gudhi/periodic_cubical_complex.pyx +++ b/src/python/gudhi/periodic_cubical_complex.pyx @@ -81,9 +81,7 @@ cdef class PeriodicCubicalComplex: periodic_dimensions=None, perseus_file=''): if ((dimensions is not None) and (top_dimensional_cells is not None) and (periodic_dimensions is not None) and (perseus_file == '')): - self.thisptr = new Periodic_cubical_complex_base_interface(dimensions, - top_dimensional_cells, - periodic_dimensions) + self._construct_from_cells(dimensions, top_dimensional_cells, periodic_dimensions) elif ((dimensions is None) and (top_dimensional_cells is not None) and (periodic_dimensions is not None) and (perseus_file == '')): top_dimensional_cells = np.array(top_dimensional_cells, @@ -91,13 +89,11 @@ cdef class PeriodicCubicalComplex: order = 'F') dimensions = top_dimensional_cells.shape top_dimensional_cells = top_dimensional_cells.ravel(order='F') - self.thisptr = new Periodic_cubical_complex_base_interface(dimensions, - top_dimensional_cells, - periodic_dimensions) + self._construct_from_cells(dimensions, top_dimensional_cells, periodic_dimensions) elif ((dimensions is None) and (top_dimensional_cells is None) and (periodic_dimensions is None) and (perseus_file != '')): if os.path.isfile(perseus_file): - self.thisptr = new Periodic_cubical_complex_base_interface(perseus_file.encode('utf-8')) + self._construct_from_file(perseus_file.encode('utf-8')) else: print("file " + perseus_file + " not found.", file=sys.stderr) else: @@ -106,6 +102,14 @@ cdef class PeriodicCubicalComplex: "top_dimensional_cells and periodic_dimensions or from " "a Perseus-style file name.", file=sys.stderr) + def _construct_from_cells(self, vector[unsigned] dimensions, vector[double] top_dimensional_cells, vector[bool] periodic_dimensions): + with nogil: + self.thisptr = new Periodic_cubical_complex_base_interface(dimensions, top_dimensional_cells, periodic_dimensions) + + def _construct_from_file(self, string filename): + with nogil: + self.thisptr = new Periodic_cubical_complex_base_interface(filename) + def __dealloc__(self): if self.thisptr != NULL: del self.thisptr -- cgit v1.2.3 From c156309dfd00c6180f2fd2dc03be159fd21c2626 Mon Sep 17 00:00:00 2001 From: Marc Glisse Date: Sun, 17 May 2020 23:32:21 +0200 Subject: One more nogil in cubical --- src/python/gudhi/cubical_complex.pyx | 3 ++- src/python/gudhi/periodic_cubical_complex.pyx | 3 ++- 2 files changed, 4 insertions(+), 2 deletions(-) (limited to 'src') diff --git a/src/python/gudhi/cubical_complex.pyx b/src/python/gudhi/cubical_complex.pyx index 0068f2ff..3ace2517 100644 --- a/src/python/gudhi/cubical_complex.pyx +++ b/src/python/gudhi/cubical_complex.pyx @@ -215,7 +215,8 @@ cdef class CubicalComplex: cdef vector[vector[int]] persistence_result output = [[],[]] - persistence_result = self.pcohptr.cofaces_of_cubical_persistence_pairs() + with nogil: + persistence_result = self.pcohptr.cofaces_of_cubical_persistence_pairs() pr = np.array(persistence_result) ess_ind = np.argwhere(pr[:,2] == -1)[:,0] diff --git a/src/python/gudhi/periodic_cubical_complex.pyx b/src/python/gudhi/periodic_cubical_complex.pyx index 11e1766c..bed55101 100644 --- a/src/python/gudhi/periodic_cubical_complex.pyx +++ b/src/python/gudhi/periodic_cubical_complex.pyx @@ -214,7 +214,8 @@ cdef class PeriodicCubicalComplex: cdef vector[vector[int]] persistence_result if self.pcohptr != NULL: output = [[],[]] - persistence_result = self.pcohptr.cofaces_of_cubical_persistence_pairs() + with nogil: + persistence_result = self.pcohptr.cofaces_of_cubical_persistence_pairs() pr = np.array(persistence_result) ess_ind = np.argwhere(pr[:,2] == -1)[:,0] -- cgit v1.2.3 From beadbbbefa1f8f30233a534b6c9cdf11ffb65f93 Mon Sep 17 00:00:00 2001 From: ROUVREAU Vincent Date: Mon, 18 May 2020 07:40:59 +0200 Subject: When Reviewing dependencies, I missed this one --- src/python/doc/wasserstein_distance_sum.inc | 20 +++++++++----------- 1 file changed, 9 insertions(+), 11 deletions(-) (limited to 'src') diff --git a/src/python/doc/wasserstein_distance_sum.inc b/src/python/doc/wasserstein_distance_sum.inc index f9308e5e..c41de017 100644 --- a/src/python/doc/wasserstein_distance_sum.inc +++ b/src/python/doc/wasserstein_distance_sum.inc @@ -1,14 +1,12 @@ .. table:: :widths: 30 40 30 - +-----------------------------------------------------------------+----------------------------------------------------------------------+------------------------------------------------------------------+ - | .. figure:: | The q-Wasserstein distance measures the similarity between two | :Author: Theo Lacombe | - | ../../doc/Bottleneck_distance/perturb_pd.png | persistence diagrams using the sum of all edges lengths (instead of | | - | :figclass: align-center | the maximum). It allows to define sophisticated objects such as | :Since: GUDHI 3.1.0 | - | | barycenters of a family of persistence diagrams. | | - | | | :License: MIT | - | | | | - | | | :Requires: Python Optimal Transport (POT) :math:`\geq` 0.5.1 | - +-----------------------------------------------------------------+----------------------------------------------------------------------+------------------------------------------------------------------+ - | * :doc:`wasserstein_distance_user` | | - +-----------------------------------------------------------------+-----------------------------------------------------------------------------------------------------------------------------------------+ + +-----------------------------------------------------------------+----------------------------------------------------------------------+-----------------------------------------+ + | .. figure:: | The q-Wasserstein distance measures the similarity between two | :Author: Theo Lacombe, Marc Glisse | + | ../../doc/Bottleneck_distance/perturb_pd.png | persistence diagrams using the sum of all edges lengths (instead of | | + | :figclass: align-center | the maximum). It allows to define sophisticated objects such as | :Since: GUDHI 3.1.0 | + | | barycenters of a family of persistence diagrams. | | + | | | :License: MIT, BSD-3-Clause | + +-----------------------------------------------------------------+----------------------------------------------------------------------+-----------------------------------------+ + | * :doc:`wasserstein_distance_user` | | + +-----------------------------------------------------------------+----------------------------------------------------------------------------------------------------------------+ -- cgit v1.2.3 From 0583c72cd729fb2d4a3e704949051e98b24726b3 Mon Sep 17 00:00:00 2001 From: MathieuCarriere Date: Wed, 20 May 2020 14:13:49 -0400 Subject: fix bug --- src/python/gudhi/cubical_complex.pyx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) (limited to 'src') diff --git a/src/python/gudhi/cubical_complex.pyx b/src/python/gudhi/cubical_complex.pyx index ca979eda..ff0a9ed8 100644 --- a/src/python/gudhi/cubical_complex.pyx +++ b/src/python/gudhi/cubical_complex.pyx @@ -209,14 +209,14 @@ cdef class CubicalComplex: ess_ind = np.argwhere(pr[:,2] == -1)[:,0] ess = pr[ess_ind] - max_h = max(ess[:,0])+1 + max_h = max(ess[:,0])+1 if len(ess) > 0 else 0 for h in range(max_h): hidxs = np.argwhere(ess[:,0] == h)[:,0] output[1].append(ess[hidxs][:,1]) reg_ind = np.setdiff1d(np.array(range(len(pr))), ess_ind) reg = pr[reg_ind] - max_h = max(reg[:,0])+1 + max_h = max(reg[:,0])+1 if len(reg) > 0 else 0 for h in range(max_h): hidxs = np.argwhere(reg[:,0] == h)[:,0] output[0].append(reg[hidxs][:,1:]) -- cgit v1.2.3 From d00bd43a7f3651e67c1572fea38550d310a223ec Mon Sep 17 00:00:00 2001 From: ROUVREAU Vincent Date: Thu, 21 May 2020 09:20:37 +0200 Subject: periodic fix and adds corresponding unitary tests for bug fix --- src/python/gudhi/periodic_cubical_complex.pyx | 42 +++++++++++++-------------- src/python/test/test_cubical_complex.py | 17 +++++++++++ 2 files changed, 37 insertions(+), 22 deletions(-) (limited to 'src') diff --git a/src/python/gudhi/periodic_cubical_complex.pyx b/src/python/gudhi/periodic_cubical_complex.pyx index bed55101..d353d2af 100644 --- a/src/python/gudhi/periodic_cubical_complex.pyx +++ b/src/python/gudhi/periodic_cubical_complex.pyx @@ -211,29 +211,27 @@ cdef class PeriodicCubicalComplex: The indices of the arrays in the list correspond to the homological dimensions, and the integers of each row in each array correspond to: (index of positive top-dimensional cell). """ + assert self.pcohptr != NULL, "compute_persistence() must be called before cofaces_of_persistence_pairs()" cdef vector[vector[int]] persistence_result - if self.pcohptr != NULL: - output = [[],[]] - with nogil: - persistence_result = self.pcohptr.cofaces_of_cubical_persistence_pairs() - pr = np.array(persistence_result) - - ess_ind = np.argwhere(pr[:,2] == -1)[:,0] - ess = pr[ess_ind] - max_h = max(ess[:,0])+1 - for h in range(max_h): - hidxs = np.argwhere(ess[:,0] == h)[:,0] - output[1].append(ess[hidxs][:,1]) - - reg_ind = np.setdiff1d(np.array(range(len(pr))), ess_ind) - reg = pr[reg_ind] - max_h = max(reg[:,0])+1 - for h in range(max_h): - hidxs = np.argwhere(reg[:,0] == h)[:,0] - output[0].append(reg[hidxs][:,1:]) - else: - print("cofaces_of_persistence_pairs function requires persistence function" - " to be launched first.") + + output = [[],[]] + with nogil: + persistence_result = self.pcohptr.cofaces_of_cubical_persistence_pairs() + pr = np.array(persistence_result) + + ess_ind = np.argwhere(pr[:,2] == -1)[:,0] + ess = pr[ess_ind] + max_h = max(ess[:,0])+1 if len(ess) > 0 else 0 + for h in range(max_h): + hidxs = np.argwhere(ess[:,0] == h)[:,0] + output[1].append(ess[hidxs][:,1]) + + reg_ind = np.setdiff1d(np.array(range(len(pr))), ess_ind) + reg = pr[reg_ind] + max_h = max(reg[:,0])+1 if len(reg) > 0 else 0 + for h in range(max_h): + hidxs = np.argwhere(reg[:,0] == h)[:,0] + output[0].append(reg[hidxs][:,1:]) return output def betti_numbers(self): diff --git a/src/python/test/test_cubical_complex.py b/src/python/test/test_cubical_complex.py index 5c59db8f..d0e4e9e8 100755 --- a/src/python/test/test_cubical_complex.py +++ b/src/python/test/test_cubical_complex.py @@ -157,3 +157,20 @@ def test_cubical_generators(): assert np.array_equal(g[0][0], np.empty(shape=[0,2])) assert np.array_equal(g[0][1], np.array([[7, 4]])) assert np.array_equal(g[1][0], np.array([8])) + +def test_cubical_cofaces_of_persistence_pairs_when_pd_has_no_paired_birth_and_death(): + cubCpx = CubicalComplex(dimensions=[1,2], top_dimensional_cells=[0.0, 1.0]) + Diag = cubCpx.persistence(homology_coeff_field=2, min_persistence=0) + pairs = cubCpx.cofaces_of_persistence_pairs() + assert pairs[0] == [] + assert np.array_equal(pairs[1][0], np.array([0])) + +def test_periodic_cofaces_of_persistence_pairs_when_pd_has_no_paired_birth_and_death(): + perCubCpx = PeriodicCubicalComplex(dimensions=[1,2], top_dimensional_cells=[0.0, 1.0], + periodic_dimensions=[True, True]) + Diag = perCubCpx.persistence(homology_coeff_field=2, min_persistence=0) + pairs = perCubCpx.cofaces_of_persistence_pairs() + assert pairs[0] == [] + assert np.array_equal(pairs[1][0], np.array([0])) + assert np.array_equal(pairs[1][1], np.array([0, 1])) + assert np.array_equal(pairs[1][2], np.array([1])) -- cgit v1.2.3