diff options
| -rw-r--r-- | src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_base.h | 21 | ||||
| -rw-r--r-- | src/python/gudhi/cubical_complex.pyx | 30 | ||||
| -rw-r--r-- | src/python/gudhi/periodic_cubical_complex.pyx | 30 | ||||
| -rw-r--r-- | src/python/include/Persistent_cohomology_interface.h | 38 | ||||
| -rwxr-xr-x | src/python/test/test_cubical_complex.py | 5 |
5 files changed, 123 insertions, 1 deletions
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 1eb77c9c..6441c129 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,6 +110,14 @@ class Bitmap_cubical_complex_base { virtual inline std::vector<std::size_t> 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 * between \f$A\f$ and \f$B\f$ is the integer with which \f$B\f$ appears in the boundary of \f$A\f$. @@ -603,6 +611,19 @@ void Bitmap_cubical_complex_base<T>::setup_bitmap_based_on_top_dimensional_cells } template <typename T> +int Bitmap_cubical_complex_base<T>::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 <typename T> Bitmap_cubical_complex_base<T>::Bitmap_cubical_complex_base(const std::vector<unsigned>& sizes_in_following_directions, const std::vector<T>& top_dimensional_cells) { this->setup_bitmap_based_on_top_dimensional_cells_list(sizes_in_following_directions, top_dimensional_cells); diff --git a/src/python/gudhi/cubical_complex.pyx b/src/python/gudhi/cubical_complex.pyx index 007abcb6..69d0f0b6 100644 --- a/src/python/gudhi/cubical_complex.pyx +++ b/src/python/gudhi/cubical_complex.pyx @@ -37,6 +37,7 @@ cdef extern from "Persistent_cohomology_interface.h" namespace "Gudhi": 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) @@ -170,6 +171,35 @@ cdef class CubicalComplex: self.compute_persistence(homology_coeff_field, min_persistence) 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 + 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/periodic_cubical_complex.pyx b/src/python/gudhi/periodic_cubical_complex.pyx index 246a3a02..78565cf8 100644 --- a/src/python/gudhi/periodic_cubical_complex.pyx +++ b/src/python/gudhi/periodic_cubical_complex.pyx @@ -34,6 +34,7 @@ cdef extern from "Persistent_cohomology_interface.h" namespace "Gudhi": 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) @@ -175,6 +176,35 @@ cdef class PeriodicCubicalComplex: self.compute_persistence(homology_coeff_field, min_persistence) 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 + 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/include/Persistent_cohomology_interface.h b/src/python/include/Persistent_cohomology_interface.h index e2b69a52..59024212 100644 --- a/src/python/include/Persistent_cohomology_interface.h +++ b/src/python/include/Persistent_cohomology_interface.h @@ -58,7 +58,6 @@ persistent_cohomology::Persistent_cohomology<FilteredComplex, persistent_cohomol cmp_intervals_by_dim_then_length cmp(stptr_); auto persistent_pairs = Base::get_persistent_pairs(); std::sort(std::begin(persistent_pairs), std::end(persistent_pairs), cmp); - std::vector<std::pair<int, std::pair<double, double>>> persistence; for (auto pair : persistent_pairs) { persistence.push_back(std::make_pair(stptr_->dimension(get<0>(pair)), @@ -68,6 +67,43 @@ persistent_cohomology::Persistent_cohomology<FilteredComplex, persistent_cohomol return persistence; } + std::vector<std::vector<int>> cofaces_of_cubical_persistence_pairs() { + + // Warning: this function is meant to be used with CubicalComplex only!! + + auto pairs = persistent_cohomology::Persistent_cohomology<FilteredComplex, + persistent_cohomology::Field_Zp>::get_persistent_pairs(); + + // Gather all top-dimensional cells and store their simplex handles + std::vector<int> 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<int, int> 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::vector<std::vector<int>> 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)); + // 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]; + } + std::vector<int> vect{ h, splx0, splx1}; + persistence_pairs.push_back(vect); + } + return persistence_pairs; + } + std::vector<std::pair<std::vector<int>, std::vector<int>>> persistence_pairs() { auto pairs = Base::get_persistent_pairs(); diff --git a/src/python/test/test_cubical_complex.py b/src/python/test/test_cubical_complex.py index fce4875c..dd7653c2 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_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]])) |