diff options
author | Marc Glisse <marc.glisse@inria.fr> | 2020-04-29 23:00:17 +0200 |
---|---|---|
committer | Marc Glisse <marc.glisse@inria.fr> | 2020-04-29 23:00:17 +0200 |
commit | 73a74011e4b5af0794f0463295beca924d32e0ee (patch) | |
tree | fb825dc77e9984594109f9ea48838d0544dfca1e /src/python/include | |
parent | 74155081bb8b3330c562d5c40d7f0a32fc188012 (diff) | |
parent | 0bba67db83f33ff608366057d9c4f005fa6a514b (diff) |
Merge remote-tracking branch 'origin/master' into dtmdensity
Diffstat (limited to 'src/python/include')
10 files changed, 161 insertions, 58 deletions
diff --git a/src/python/include/Alpha_complex_interface.h b/src/python/include/Alpha_complex_interface.h index 8614eee3..40de88f3 100644 --- a/src/python/include/Alpha_complex_interface.h +++ b/src/python/include/Alpha_complex_interface.h @@ -58,7 +58,6 @@ class Alpha_complex_interface { void create_simplex_tree(Simplex_tree_interface<>* simplex_tree, double max_alpha_square) { alpha_complex_->create_complex(*simplex_tree, max_alpha_square); - simplex_tree->initialize_filtration(); } private: diff --git a/src/python/include/Euclidean_strong_witness_complex_interface.h b/src/python/include/Euclidean_strong_witness_complex_interface.h index c1c72737..f94c51ef 100644 --- a/src/python/include/Euclidean_strong_witness_complex_interface.h +++ b/src/python/include/Euclidean_strong_witness_complex_interface.h @@ -50,12 +50,10 @@ class Euclidean_strong_witness_complex_interface { void create_simplex_tree(Gudhi::Simplex_tree<>* simplex_tree, double max_alpha_square, std::size_t limit_dimension) { witness_complex_->create_complex(*simplex_tree, max_alpha_square, limit_dimension); - simplex_tree->initialize_filtration(); } void create_simplex_tree(Gudhi::Simplex_tree<>* simplex_tree, double max_alpha_square) { witness_complex_->create_complex(*simplex_tree, max_alpha_square); - simplex_tree->initialize_filtration(); } std::vector<double> get_point(unsigned vh) { diff --git a/src/python/include/Euclidean_witness_complex_interface.h b/src/python/include/Euclidean_witness_complex_interface.h index 5d7dbdc2..4411ae79 100644 --- a/src/python/include/Euclidean_witness_complex_interface.h +++ b/src/python/include/Euclidean_witness_complex_interface.h @@ -49,12 +49,10 @@ class Euclidean_witness_complex_interface { void create_simplex_tree(Gudhi::Simplex_tree<>* simplex_tree, double max_alpha_square, std::size_t limit_dimension) { witness_complex_->create_complex(*simplex_tree, max_alpha_square, limit_dimension); - simplex_tree->initialize_filtration(); } void create_simplex_tree(Gudhi::Simplex_tree<>* simplex_tree, double max_alpha_square) { witness_complex_->create_complex(*simplex_tree, max_alpha_square); - simplex_tree->initialize_filtration(); } std::vector<double> get_point(unsigned vh) { diff --git a/src/python/include/Nerve_gic_interface.h b/src/python/include/Nerve_gic_interface.h index 5e7f8ae6..ab14c318 100644 --- a/src/python/include/Nerve_gic_interface.h +++ b/src/python/include/Nerve_gic_interface.h @@ -29,7 +29,6 @@ class Nerve_gic_interface : public Cover_complex<std::vector<double>> { public: void create_simplex_tree(Simplex_tree_interface<>* simplex_tree) { create_complex(*simplex_tree); - simplex_tree->initialize_filtration(); } void set_cover_from_Euclidean_Voronoi(int m) { set_cover_from_Voronoi(Gudhi::Euclidean_distance(), m); diff --git a/src/python/include/Persistent_cohomology_interface.h b/src/python/include/Persistent_cohomology_interface.h index 8c79e6f3..0de9bd5c 100644 --- a/src/python/include/Persistent_cohomology_interface.h +++ b/src/python/include/Persistent_cohomology_interface.h @@ -23,82 +23,162 @@ 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), - stptr_(stptr) { } - - Persistent_cohomology_interface(FilteredComplex* stptr, bool persistence_dim_max) - : persistent_cohomology::Persistent_cohomology<FilteredComplex, - persistent_cohomology::Field_Zp>(*stptr, persistence_dim_max), + Persistent_cohomology_interface(FilteredComplex* stptr, bool persistence_dim_max=false) + : 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); + // TODO: move to the constructors? + void compute_persistence(int homology_coeff_field, double min_persistence) { + Base::init_coefficients(homology_coeff_field); + Base::compute_persistent_cohomology(min_persistence); + } + std::vector<std::pair<int, std::pair<double, double>>> get_persistence() { std::vector<std::pair<int, std::pair<double, double>>> persistence; + auto const& persistent_pairs = Base::get_persistent_pairs(); + 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() { + 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); + 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() { + 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(); + 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_; diff --git a/src/python/include/Rips_complex_interface.h b/src/python/include/Rips_complex_interface.h index a66b0e5b..d98b0226 100644 --- a/src/python/include/Rips_complex_interface.h +++ b/src/python/include/Rips_complex_interface.h @@ -53,7 +53,6 @@ class Rips_complex_interface { rips_complex_->create_complex(*simplex_tree, dim_max); else sparse_rips_complex_->create_complex(*simplex_tree, dim_max); - simplex_tree->initialize_filtration(); } private: diff --git a/src/python/include/Simplex_tree_interface.h b/src/python/include/Simplex_tree_interface.h index 4a7062d6..56d7c41d 100644 --- a/src/python/include/Simplex_tree_interface.h +++ b/src/python/include/Simplex_tree_interface.h @@ -16,8 +16,6 @@ #include <gudhi/Simplex_tree.h> #include <gudhi/Points_off_io.h> -#include "Persistent_cohomology_interface.h" - #include <iostream> #include <vector> #include <utility> // std::pair @@ -37,18 +35,25 @@ class Simplex_tree_interface : public Simplex_tree<SimplexTreeOptions> { using Filtered_simplices = std::vector<Simplex_and_filtration>; using Skeleton_simplex_iterator = typename Base::Skeleton_simplex_iterator; using Complex_simplex_iterator = typename Base::Complex_simplex_iterator; + using Extended_filtration_data = typename Base::Extended_filtration_data; public: - bool find_simplex(const Simplex& vh) { - return (Base::find(vh) != Base::null_simplex()); + + Extended_filtration_data efd; + + bool find_simplex(const Simplex& simplex) { + return (Base::find(simplex) != Base::null_simplex()); } - void assign_simplex_filtration(const Simplex& vh, Filtration_value filtration) { - Base::assign_filtration(Base::find(vh), filtration); + void assign_simplex_filtration(const Simplex& simplex, Filtration_value filtration) { + Base::assign_filtration(Base::find(simplex), filtration); + Base::clear_filtration(); } bool insert(const Simplex& simplex, Filtration_value filtration = 0) { Insertion_result result = Base::insert_simplex_and_subfaces(simplex, filtration); + if (result.first != Base::null_simplex()) + Base::clear_filtration(); return (result.second); } @@ -82,7 +87,7 @@ class Simplex_tree_interface : public Simplex_tree<SimplexTreeOptions> { void remove_maximal_simplex(const Simplex& simplex) { Base::remove_maximal_simplex(Base::find(simplex)); - Base::initialize_filtration(); + Base::clear_filtration(); } Simplex_and_filtration get_simplex_and_filtration(Simplex_handle f_simplex) { @@ -117,9 +122,39 @@ class Simplex_tree_interface : public Simplex_tree<SimplexTreeOptions> { return cofaces; } - void create_persistence(Gudhi::Persistent_cohomology_interface<Base>* pcoh) { - Base::initialize_filtration(); - pcoh = new Gudhi::Persistent_cohomology_interface<Base>(*this); + void compute_extended_filtration() { + this->efd = this->extend_filtration(); + return; + } + + std::vector<std::vector<std::pair<int, std::pair<Filtration_value, Filtration_value>>>> compute_extended_persistence_subdiagrams(const std::vector<std::pair<int, std::pair<Filtration_value, Filtration_value>>>& dgm, Filtration_value min_persistence){ + std::vector<std::vector<std::pair<int, std::pair<Filtration_value, Filtration_value>>>> new_dgm(4); + for (unsigned int i = 0; i < dgm.size(); i++){ + std::pair<Filtration_value, Extended_simplex_type> px = this->decode_extended_filtration(dgm[i].second.first, this->efd); + std::pair<Filtration_value, Extended_simplex_type> py = this->decode_extended_filtration(dgm[i].second.second, this->efd); + std::pair<int, std::pair<Filtration_value, Filtration_value>> pd_point = std::make_pair(dgm[i].first, std::make_pair(px.first, py.first)); + if(std::abs(px.first - py.first) > min_persistence){ + //Ordinary + if (px.second == Extended_simplex_type::UP && py.second == Extended_simplex_type::UP){ + new_dgm[0].push_back(pd_point); + } + // Relative + else if (px.second == Extended_simplex_type::DOWN && py.second == Extended_simplex_type::DOWN){ + new_dgm[1].push_back(pd_point); + } + else{ + // Extended+ + if (px.first < py.first){ + new_dgm[2].push_back(pd_point); + } + //Extended- + else{ + new_dgm[3].push_back(pd_point); + } + } + } + } + return new_dgm; } // Iterator over the simplex tree diff --git a/src/python/include/Strong_witness_complex_interface.h b/src/python/include/Strong_witness_complex_interface.h index cda5b514..e9ab0c7b 100644 --- a/src/python/include/Strong_witness_complex_interface.h +++ b/src/python/include/Strong_witness_complex_interface.h @@ -41,13 +41,11 @@ class Strong_witness_complex_interface { void create_simplex_tree(Simplex_tree_interface<>* simplex_tree, double max_alpha_square, std::size_t limit_dimension) { witness_complex_->create_complex(*simplex_tree, max_alpha_square, limit_dimension); - simplex_tree->initialize_filtration(); } void create_simplex_tree(Simplex_tree_interface<>* simplex_tree, double max_alpha_square) { witness_complex_->create_complex(*simplex_tree, max_alpha_square); - simplex_tree->initialize_filtration(); } private: diff --git a/src/python/include/Tangential_complex_interface.h b/src/python/include/Tangential_complex_interface.h index 698226cc..b1afce94 100644 --- a/src/python/include/Tangential_complex_interface.h +++ b/src/python/include/Tangential_complex_interface.h @@ -90,7 +90,6 @@ class Tangential_complex_interface { void create_simplex_tree(Simplex_tree<>* simplex_tree) { tangential_complex_->create_complex<Gudhi::Simplex_tree<Gudhi::Simplex_tree_options_full_featured>>(*simplex_tree); - simplex_tree->initialize_filtration(); } void set_max_squared_edge_length(double max_squared_edge_length) { diff --git a/src/python/include/Witness_complex_interface.h b/src/python/include/Witness_complex_interface.h index 45e14253..76947e53 100644 --- a/src/python/include/Witness_complex_interface.h +++ b/src/python/include/Witness_complex_interface.h @@ -41,13 +41,11 @@ class Witness_complex_interface { void create_simplex_tree(Simplex_tree_interface<>* simplex_tree, double max_alpha_square, std::size_t limit_dimension) { witness_complex_->create_complex(*simplex_tree, max_alpha_square, limit_dimension); - simplex_tree->initialize_filtration(); } void create_simplex_tree(Simplex_tree_interface<>* simplex_tree, double max_alpha_square) { witness_complex_->create_complex(*simplex_tree, max_alpha_square); - simplex_tree->initialize_filtration(); } private: |