diff options
Diffstat (limited to 'src/python/include/Simplex_tree_interface.h')
| -rw-r--r-- | src/python/include/Simplex_tree_interface.h | 97 |
1 files changed, 71 insertions, 26 deletions
diff --git a/src/python/include/Simplex_tree_interface.h b/src/python/include/Simplex_tree_interface.h index 56d7c41d..0317ea39 100644 --- a/src/python/include/Simplex_tree_interface.h +++ b/src/python/include/Simplex_tree_interface.h @@ -15,10 +15,13 @@ #include <gudhi/distance_functions.h> #include <gudhi/Simplex_tree.h> #include <gudhi/Points_off_io.h> +#include <gudhi/Flag_complex_edge_collapser.h> #include <iostream> #include <vector> #include <utility> // std::pair +#include <tuple> +#include <iterator> // for std::distance namespace Gudhi { @@ -36,6 +39,10 @@ class Simplex_tree_interface : public Simplex_tree<SimplexTreeOptions> { 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; + using Boundary_simplex_iterator = typename Base::Boundary_simplex_iterator; + using Siblings = typename Base::Siblings; + using Node = typename Base::Node; + typedef bool (*blocker_func_t)(Simplex simplex, void *user_data); public: @@ -57,6 +64,30 @@ class Simplex_tree_interface : public Simplex_tree<SimplexTreeOptions> { return (result.second); } + void insert_matrix(double* filtrations, int n, int stride0, int stride1, double max_filtration) { + // We could delegate to insert_graph, but wrapping the matrix in a graph interface is too much work, + // and this is a bit more efficient. + auto& rm = this->root()->members_; + for(int i=0; i<n; ++i) { + char* p = reinterpret_cast<char*>(filtrations) + i * stride0; + double fv = *reinterpret_cast<double*>(p + i * stride1); + if(fv > max_filtration) continue; + auto sh = rm.emplace_hint(rm.end(), i, Node(this->root(), fv)); + Siblings* children = nullptr; + // Should we make a first pass to count the number of edges so we can reserve the right space? + for(int j=i+1; j<n; ++j) { + double fe = *reinterpret_cast<double*>(p + j * stride1); + if(fe > max_filtration) continue; + if(!children) { + children = new Siblings(this->root(), i); + sh->second.assign_children(children); + } + children->members().emplace_hint(children->members().end(), j, Node(children, fe)); + } + } + + } + // Do not interface this function, only used in alpha complex interface for complex creation bool insert_simplex(const Simplex& simplex, Filtration_value filtration = 0) { Insertion_result result = Base::insert_simplex(simplex, filtration); @@ -127,34 +158,39 @@ class Simplex_tree_interface : public Simplex_tree<SimplexTreeOptions> { 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); - } - } + Simplex_tree_interface* collapse_edges(int nb_collapse_iteration) { + using Filtered_edge = std::tuple<Vertex_handle, Vertex_handle, Filtration_value>; + std::vector<Filtered_edge> edges; + for (Simplex_handle sh : Base::skeleton_simplex_range(1)) { + if (Base::dimension(sh) == 1) { + typename Base::Simplex_vertex_range rg = Base::simplex_vertex_range(sh); + auto vit = rg.begin(); + Vertex_handle v = *vit; + Vertex_handle w = *++vit; + edges.emplace_back(v, w, Base::filtration(sh)); } } - return new_dgm; + + for (int iteration = 0; iteration < nb_collapse_iteration; iteration++) { + edges = Gudhi::collapse::flag_complex_collapse_edges(std::move(edges)); + } + Simplex_tree_interface* collapsed_stree_ptr = new Simplex_tree_interface(); + // Copy the original 0-skeleton + for (Simplex_handle sh : Base::skeleton_simplex_range(0)) { + collapsed_stree_ptr->insert({*(Base::simplex_vertex_range(sh).begin())}, Base::filtration(sh)); + } + // Insert remaining edges + for (auto remaining_edge : edges) { + collapsed_stree_ptr->insert({std::get<0>(remaining_edge), std::get<1>(remaining_edge)}, std::get<2>(remaining_edge)); + } + return collapsed_stree_ptr; + } + + void expansion_with_blockers_callback(int dimension, blocker_func_t user_func, void *user_data) { + Base::expansion_with_blockers(dimension, [&](Simplex_handle sh){ + Simplex simplex(Base::simplex_vertex_range(sh).begin(), Base::simplex_vertex_range(sh).end()); + return user_func(simplex, user_data); + }); } // Iterator over the simplex tree @@ -188,6 +224,15 @@ class Simplex_tree_interface : public Simplex_tree<SimplexTreeOptions> { // this specific case works because the range is just a pair of iterators - won't work if range was a vector return Base::skeleton_simplex_range(dimension).end(); } + + std::pair<Boundary_simplex_iterator, Boundary_simplex_iterator> get_boundary_iterators(const Simplex& simplex) { + auto bd_sh = Base::find(simplex); + if (bd_sh == Base::null_simplex()) + throw std::runtime_error("simplex not found - cannot find boundaries"); + // this specific case works because the range is just a pair of iterators - won't work if range was a vector + auto boundary_srange = Base::boundary_simplex_range(bd_sh); + return std::make_pair(boundary_srange.begin(), boundary_srange.end()); + } }; } // namespace Gudhi |