diff options
Diffstat (limited to 'src/python')
-rw-r--r-- | src/python/gudhi/simplex_tree.pxd | 2 | ||||
-rw-r--r-- | src/python/gudhi/simplex_tree.pyx | 64 | ||||
-rw-r--r-- | src/python/include/Persistent_cohomology_interface.h | 146 |
3 files changed, 181 insertions, 31 deletions
diff --git a/src/python/gudhi/simplex_tree.pxd b/src/python/gudhi/simplex_tree.pxd index 82f155de..44789365 100644 --- a/src/python/gudhi/simplex_tree.pxd +++ b/src/python/gudhi/simplex_tree.pxd @@ -75,3 +75,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(double) + pair[vector[vector[int]], vector[vector[int]]] flag_generators(double) diff --git a/src/python/gudhi/simplex_tree.pyx b/src/python/gudhi/simplex_tree.pyx index c01cc905..faa9f9d8 100644 --- a/src/python/gudhi/simplex_tree.pyx +++ b/src/python/gudhi/simplex_tree.pyx @@ -9,6 +9,7 @@ from cython.operator import dereference, preincrement from libc.stdint cimport intptr_t +import numpy from numpy import array as np_array cimport simplex_tree @@ -405,7 +406,7 @@ cdef class SimplexTree: :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. + Set 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 @@ -524,3 +525,64 @@ cdef class SimplexTree: else: print("intervals_in_dim function requires persistence function" " to be launched first.") + + def lower_star_persistence_generators(self, min_persistence=0.): + """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. + + :param min_persistence: The minimum persistence value to take into + account (strictly greater than min_persistence). Default value is + 0.0. + Set min_persistence to -1.0 to see all values. + :type min_persistence: float. + :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: lower_star_persistence_generators requires that `persistence()` be called first. + """ + if self.pcohptr != NULL: + gen = self.pcohptr.lower_star_generators(min_persistence) + 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, min_persistence=0.): + """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. + + :param min_persistence: The minimum persistence value to take into + account (strictly greater than min_persistence). Default value is + 0.0. + Set min_persistence to -1.0 to see all values. + :type min_persistence: float. + :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[numpy.array[int] of shape (n,3), List[numpy.array[int] of shape (m,4)], numpy.array[int] of shape (l,), List[numpy.array[int] of shape (k,2)]] + + :note: flag_persistence_generators requires that `persistence()` be called first. + """ + if self.pcohptr != NULL: + gen = self.pcohptr.flag_generators(min_persistence) + if len(gen.first) == 0: + normal0 = numpy.empty((0,3)) + normals = [] + 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 = numpy.empty(0) + infinites = [] + else: + l = iter(gen.second) + infinite0 = np_array(next(l)) + infinites = [np_array(d).reshape(-1,2) for d in l] + + return (normal0, normals, infinite0, infinites) + else: + print("flag_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..3ce40af5 100644 --- a/src/python/include/Persistent_cohomology_interface.h +++ b/src/python/include/Persistent_cohomology_interface.h @@ -23,82 +23,168 @@ 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()); + std::vector<int> birth; + std::vector<int> death; for (auto pair : pairs) { - std::vector<int> birth; + birth.clear(); if (get<0>(pair) != stptr_->null_simplex()) { for (auto vertex : stptr_->simplex_vertex_range(get<0>(pair))) { birth.push_back(vertex); } } - std::vector<int> death; + death.clear(); 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(birth, death); } return persistence_pairs; } + // TODO: (possibly at the python level) + // - 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(double min_persistence) { + 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); + if(stptr_->filtration(t) - stptr_->filtration(s) <= min_persistence) + continue; + 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); + auto& d = diags[dim]; + d.insert(d.end(), { v, w }); + } + } + return out; + } + + // An alternative, to avoid those different sizes, would be to "pad" vertex generator v as (v, v) or (v, -1). When using it as index, this corresponds to adding the vertex filtration values either on the diagonal of the distance matrix, or as an extra row or column. + // We could also merge the vectors for different dimensions into a single one, with an extra column for the dimension (converted to type double). + Generators flag_generators(double min_persistence) { + 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); + if(stptr_->filtration(t) - stptr_->filtration(s) <= min_persistence) + continue; + 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(); + auto& d = diagsinf[dim]; + d.insert(d.end(), { v1, 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(); + auto& d = diags[0]; + d.insert(d.end(), { v, w1, 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(); + auto& d = diags[dim]; + d.insert(d.end(), { v1, v2, w1, w2 }); + } + } + } + return out; + } + private: // A copy FilteredComplex* stptr_; |