diff options
author | Marc Glisse <marc.glisse@inria.fr> | 2020-03-07 01:08:10 +0100 |
---|---|---|
committer | Marc Glisse <marc.glisse@inria.fr> | 2020-03-07 01:08:10 +0100 |
commit | 64199fd8037556f135f90102ba8270cccf9d3e60 (patch) | |
tree | eba5ee334253e2afc7291330f755283f592e911d | |
parent | 8bd39f74f69e8fcb662873e0e045c953b814f28f (diff) |
persistence generators for lower-star and flag filtrations
-rw-r--r-- | src/python/gudhi/simplex_tree.pxd | 2 | ||||
-rw-r--r-- | src/python/gudhi/simplex_tree.pyx | 55 | ||||
-rw-r--r-- | src/python/include/Persistent_cohomology_interface.h | 138 |
3 files changed, 167 insertions, 28 deletions
diff --git a/src/python/gudhi/simplex_tree.pxd b/src/python/gudhi/simplex_tree.pxd index 96d14079..4e435c67 100644 --- a/src/python/gudhi/simplex_tree.pxd +++ b/src/python/gudhi/simplex_tree.pxd @@ -53,3 +53,5 @@ cdef extern from "Persistent_cohomology_interface.h" namespace "Gudhi": vector[pair[double,double]] intervals_in_dimension(int dimension) void write_output_diagram(string diagram_file_name) vector[pair[vector[int], vector[int]]] persistence_pairs() + pair[vector[vector[int]], vector[vector[int]]] lower_star_generators() + pair[vector[vector[int]], vector[vector[int]]] flag_generators() diff --git a/src/python/gudhi/simplex_tree.pyx b/src/python/gudhi/simplex_tree.pyx index b18627c4..1c9b9cf1 100644 --- a/src/python/gudhi/simplex_tree.pyx +++ b/src/python/gudhi/simplex_tree.pyx @@ -514,3 +514,58 @@ cdef class SimplexTree: else: print("intervals_in_dim function requires persistence function" " to be launched first.") + + def lower_star_persistence_generators(self): + """Assuming this is a lower-star filtration, this function returns the persistence pairs, + where each simplex is replaced with the vertex that gave it its filtration value. + + :returns: first the regular persistence pairs, grouped by dimension, with one vertex per extremity, + and second the essential features, grouped by dimension, with one vertex each + :rtype: Tuple[List[numpy.array[int] of shape (n,2)], List[numpy.array[int] of shape (m,)]] + + :note: intervals_in_dim function requires + :func:`persistence()<gudhi.SimplexTree.persistence>` + function to be launched first. + """ + if self.pcohptr != NULL: + gen = self.pcohptr.lower_star_generators() + normal = [np_array(d).reshape(-1,2) for d in gen.first] + infinite = [np_array(d) for d in gen.second] + return (normal, infinite) + else: + print("lower_star_persistence_generators() requires that persistence() be called first.") + + def flag_persistence_generators(self): + """Assuming this is a flag complex, this function returns the persistence pairs, + where each simplex is replaced with the vertices of the edges that gave it its filtration value. + + :returns: first the regular persistence pairs of dimension 0, with one vertex for birth and two for death; + then the other regular persistence pairs, grouped by dimension, with 2 vertices per extremity; + then the connected components, with one vertex each; + finally the other essential features, grouped by dimension, with 2 vertices for birth. + :rtype: Tuple[List[numpy.array[int] of shape (n,3)], List[numpy.array[int] of shape (m,4)], List[numpy.array[int] of shape (l,)], List[numpy.array[int] of shape (k,2)]] + + :note: intervals_in_dim function requires + :func:`persistence()<gudhi.SimplexTree.persistence>` + function to be launched first. + """ + if self.pcohptr != NULL: + gen = self.pcohptr.flag_generators() + if len(gen.first) == 0: + normal0 = np_array([]) + normals = np_array([]) + else: + l = iter(gen.first) + normal0 = np_array(next(l)).reshape(-1,3) + normals = [np_array(d).reshape(-1,4) for d in l] + if len(gen.second) == 0: + infinite0 = np_array([]) + infinites = np_array([]) + else: + l = iter(gen.second) + infinite0 = np_array(next(l)) + infinites = [np_array(d).reshape(-1,3) for d in l] + + return (normal0, normals, infinite0, infinites) + else: + print("lower_star_persistence_generators() requires that persistence() be called first.") diff --git a/src/python/include/Persistent_cohomology_interface.h b/src/python/include/Persistent_cohomology_interface.h index 8c79e6f3..6e9aac52 100644 --- a/src/python/include/Persistent_cohomology_interface.h +++ b/src/python/include/Persistent_cohomology_interface.h @@ -23,61 +23,55 @@ template<class FilteredComplex> class Persistent_cohomology_interface : public persistent_cohomology::Persistent_cohomology<FilteredComplex, persistent_cohomology::Field_Zp> { private: + typedef persistent_cohomology::Persistent_cohomology<FilteredComplex, persistent_cohomology::Field_Zp> Base; /* * Compare two intervals by dimension, then by length. */ struct cmp_intervals_by_dim_then_length { - explicit cmp_intervals_by_dim_then_length(FilteredComplex * sc) - : sc_(sc) { } - template<typename Persistent_interval> bool operator()(const Persistent_interval & p1, const Persistent_interval & p2) { - if (sc_->dimension(get < 0 > (p1)) == sc_->dimension(get < 0 > (p2))) - return (sc_->filtration(get < 1 > (p1)) - sc_->filtration(get < 0 > (p1)) - > sc_->filtration(get < 1 > (p2)) - sc_->filtration(get < 0 > (p2))); + if (std::get<0>(p1) == std::get<0>(p2)) { + auto& i1 = std::get<1>(p1); + auto& i2 = std::get<1>(p2); + return std::get<1>(i1) - std::get<0>(i1) > std::get<1>(i2) - std::get<0>(i2); + } else - return (sc_->dimension(get < 0 > (p1)) > sc_->dimension(get < 0 > (p2))); + return (std::get<0>(p1) > std::get<0>(p2)); + // Why does this sort by decreasing dimension? } - FilteredComplex* sc_; }; public: Persistent_cohomology_interface(FilteredComplex* stptr) - : persistent_cohomology::Persistent_cohomology<FilteredComplex, persistent_cohomology::Field_Zp>(*stptr), + : Base(*stptr), stptr_(stptr) { } Persistent_cohomology_interface(FilteredComplex* stptr, bool persistence_dim_max) - : persistent_cohomology::Persistent_cohomology<FilteredComplex, - persistent_cohomology::Field_Zp>(*stptr, persistence_dim_max), + : Base(*stptr, persistence_dim_max), stptr_(stptr) { } std::vector<std::pair<int, std::pair<double, double>>> get_persistence(int homology_coeff_field, double min_persistence) { - persistent_cohomology::Persistent_cohomology<FilteredComplex, - persistent_cohomology::Field_Zp>::init_coefficients(homology_coeff_field); - persistent_cohomology::Persistent_cohomology<FilteredComplex, - persistent_cohomology::Field_Zp>::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<FilteredComplex, - persistent_cohomology::Field_Zp>::get_persistent_pairs(); - std::sort(std::begin(persistent_pairs), std::end(persistent_pairs), cmp); + Base::init_coefficients(homology_coeff_field); + Base::compute_persistent_cohomology(min_persistence); + auto const& persistent_pairs = Base::get_persistent_pairs(); std::vector<std::pair<int, std::pair<double, double>>> persistence; + persistence.reserve(persistent_pairs.size()); for (auto pair : persistent_pairs) { - persistence.push_back(std::make_pair(stptr_->dimension(get<0>(pair)), - std::make_pair(stptr_->filtration(get<0>(pair)), - stptr_->filtration(get<1>(pair))))); + persistence.emplace_back(stptr_->dimension(get<0>(pair)), + std::make_pair(stptr_->filtration(get<0>(pair)), + stptr_->filtration(get<1>(pair)))); } + // Custom sort and output persistence + cmp_intervals_by_dim_then_length cmp; + std::sort(std::begin(persistence), std::end(persistence), cmp); return persistence; } std::vector<std::pair<std::vector<int>, std::vector<int>>> persistence_pairs() { - auto pairs = persistent_cohomology::Persistent_cohomology<FilteredComplex, - persistent_cohomology::Field_Zp>::get_persistent_pairs(); - std::vector<std::pair<std::vector<int>, std::vector<int>>> persistence_pairs; + auto const& pairs = Base::get_persistent_pairs(); persistence_pairs.reserve(pairs.size()); for (auto pair : pairs) { std::vector<int> birth; @@ -89,16 +83,104 @@ persistent_cohomology::Persistent_cohomology<FilteredComplex, persistent_cohomol std::vector<int> death; if (get<1>(pair) != stptr_->null_simplex()) { + death.reserve(birth.size()+1); for (auto vertex : stptr_->simplex_vertex_range(get<1>(pair))) { death.push_back(vertex); } } - persistence_pairs.push_back(std::make_pair(birth, death)); + persistence_pairs.emplace_back(std::move(birth), std::move(death)); } return persistence_pairs; } + // TODO: (possibly at the python level) + // - an option to ignore intervals of length 0? + // - an option to return only some of those vectors? + typedef std::pair<std::vector<std::vector<int>>, std::vector<std::vector<int>>> Generators; + + Generators lower_star_generators() { + Generators out; + // diags[i] should be interpreted as vector<array<int,2>> + auto& diags = out.first; + // diagsinf[i] should be interpreted as vector<int> + auto& diagsinf = out.second; + for (auto pair : Base::get_persistent_pairs()) { + auto s = std::get<0>(pair); + auto t = std::get<1>(pair); + int dim = stptr_->dimension(s); + auto v = stptr_->vertex_with_same_filtration(s); + if(t == stptr_->null_simplex()) { + while(diagsinf.size() < dim+1) diagsinf.emplace_back(); + diagsinf[dim].push_back(v); + } else { + while(diags.size() < dim+1) diags.emplace_back(); + auto w = stptr_->vertex_with_same_filtration(t); + diags[dim].push_back(v); + diags[dim].push_back(w); + } + } + return out; + } + + Generators flag_generators() { + Generators out; + // diags[0] should be interpreted as vector<array<int,3>> and other diags[i] as vector<array<int,4>> + auto& diags = out.first; + // diagsinf[0] should be interpreted as vector<int> and other diagsinf[i] as vector<array<int,2>> + auto& diagsinf = out.second; + for (auto pair : Base::get_persistent_pairs()) { + auto s = std::get<0>(pair); + auto t = std::get<1>(pair); + int dim = stptr_->dimension(s); + bool infinite = t == stptr_->null_simplex(); + if(infinite) { + if(dim == 0) { + auto v = *std::begin(stptr_->simplex_vertex_range(s)); + if(diagsinf.size()==0)diagsinf.emplace_back(); + diagsinf[0].push_back(v); + } else { + auto e = stptr_->edge_with_same_filtration(s); + auto&& e_vertices = stptr_->simplex_vertex_range(e); + auto i = std::begin(e_vertices); + auto v1 = *i; + auto v2 = *++i; + GUDHI_CHECK(++i==std::end(e_vertices), "must be an edge"); + while(diagsinf.size() < dim+1) diagsinf.emplace_back(); + diagsinf[dim].push_back(v1); + diagsinf[dim].push_back(v2); + } + } else { + auto et = stptr_->edge_with_same_filtration(t); + auto&& et_vertices = stptr_->simplex_vertex_range(et); + auto it = std::begin(et_vertices); + auto w1 = *it; + auto w2 = *++it; + GUDHI_CHECK(++it==std::end(et_vertices), "must be an edge"); + if(dim == 0) { + auto v = *std::begin(stptr_->simplex_vertex_range(s)); + if(diags.size()==0)diags.emplace_back(); + diags[0].push_back(v); + diags[0].push_back(w1); + diags[0].push_back(w2); + } else { + auto es = stptr_->edge_with_same_filtration(s); + auto&& es_vertices = stptr_->simplex_vertex_range(es); + auto is = std::begin(es_vertices); + auto v1 = *is; + auto v2 = *++is; + GUDHI_CHECK(++is==std::end(es_vertices), "must be an edge"); + while(diags.size() < dim+1) diags.emplace_back(); + diags[dim].push_back(v1); + diags[dim].push_back(v2); + diags[dim].push_back(w1); + diags[dim].push_back(w2); + } + } + } + return out; + } + private: // A copy FilteredComplex* stptr_; |