diff options
author | Vincent Rouvreau <10407034+VincentRouvreau@users.noreply.github.com> | 2020-05-15 14:33:14 +0200 |
---|---|---|
committer | GitHub <noreply@github.com> | 2020-05-15 14:33:14 +0200 |
commit | 9fb4015f9ccd394146bc436d7011d7855d919837 (patch) | |
tree | bca7f533a3d9ca7bcdd4a5a5557e284f8bdd8116 /src/python | |
parent | 1efd71c502bacce375e1950e10a8112208acd0cf (diff) | |
parent | b0ae08e93fdba8a1faec56c2230b6f542653c49e (diff) |
Merge pull request #204 from MathieuCarriere/generators
Generators
Diffstat (limited to 'src/python')
-rw-r--r-- | src/python/gudhi/cubical_complex.pyx | 53 | ||||
-rw-r--r-- | src/python/gudhi/periodic_cubical_complex.pyx | 53 | ||||
-rw-r--r-- | src/python/include/Persistent_cohomology_interface.h | 44 | ||||
-rwxr-xr-x | src/python/test/test_cubical_complex.py | 10 |
4 files changed, 160 insertions, 0 deletions
diff --git a/src/python/gudhi/cubical_complex.pyx b/src/python/gudhi/cubical_complex.pyx index 007abcb6..ca979eda 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,58 @@ 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, 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). + """ + + assert self.pcohptr != NULL, "compute_persistence() must be called before cofaces_of_persistence_pairs()" + + cdef vector[vector[int]] persistence_result + 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): """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..06309772 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,58 @@ 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, 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 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 0de9bd5c..e5a3dfba 100644 --- a/src/python/include/Persistent_cohomology_interface.h +++ b/src/python/include/Persistent_cohomology_interface.h @@ -13,9 +13,11 @@ #include <gudhi/Persistent_cohomology.h> +#include <cstdlib> #include <vector> #include <utility> // for std::pair #include <algorithm> // for sort +#include <unordered_map> namespace Gudhi { @@ -67,6 +69,48 @@ persistent_cohomology::Persistent_cohomology<FilteredComplex, persistent_cohomol return persistence; } + // This function computes the top-dimensional cofaces associated to the positive and negative + // simplices of a cubical complex. The output format is a vector of vectors of three integers, + // which are [homological dimension, index of top-dimensional coface of positive simplex, + // index of top-dimensional coface of negative simplex]. If the topological feature is essential, + // then the index of top-dimensional coface of negative simplex is arbitrarily set to -1. + 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<std::size_t> 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<std::size_t, int> 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<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 + 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 (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]; + } + persistence_pairs.push_back({ h, splx0, splx1 }); + } + return persistence_pairs; + } + std::vector<std::pair<std::vector<int>, std::vector<int>>> persistence_pairs() { std::vector<std::pair<std::vector<int>, std::vector<int>>> persistence_pairs; auto const& 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..5c59db8f 100755 --- a/src/python/test/test_cubical_complex.py +++ b/src/python/test/test_cubical_complex.py @@ -147,3 +147,13 @@ 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() + 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])) |