summaryrefslogtreecommitdiff
path: root/src/python/include
diff options
context:
space:
mode:
Diffstat (limited to 'src/python/include')
-rw-r--r--src/python/include/Alpha_complex_interface.h1
-rw-r--r--src/python/include/Euclidean_strong_witness_complex_interface.h2
-rw-r--r--src/python/include/Euclidean_witness_complex_interface.h2
-rw-r--r--src/python/include/Nerve_gic_interface.h1
-rw-r--r--src/python/include/Persistent_cohomology_interface.h194
-rw-r--r--src/python/include/Rips_complex_interface.h1
-rw-r--r--src/python/include/Simplex_tree_interface.h119
-rw-r--r--src/python/include/Strong_witness_complex_interface.h2
-rw-r--r--src/python/include/Tangential_complex_interface.h1
-rw-r--r--src/python/include/Witness_complex_interface.h2
-rw-r--r--src/python/include/pybind11_diagram_utils.h39
11 files changed, 284 insertions, 80 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..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 {
@@ -23,82 +25,204 @@ 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,
+ // 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();
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 06f31341..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
@@ -33,19 +31,29 @@ class Simplex_tree_interface : public Simplex_tree<SimplexTreeOptions> {
using Simplex_handle = typename Base::Simplex_handle;
using Insertion_result = typename std::pair<Simplex_handle, bool>;
using Simplex = std::vector<Vertex_handle>;
- using Filtered_simplices = std::vector<std::pair<Simplex, Filtration_value>>;
+ using Simplex_and_filtration = std::pair<Simplex, Filtration_value>;
+ 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);
}
@@ -79,32 +87,15 @@ 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();
}
- Filtered_simplices get_filtration() {
- Base::initialize_filtration();
- Filtered_simplices filtrations;
- for (auto f_simplex : Base::filtration_simplex_range()) {
- Simplex simplex;
- for (auto vertex : Base::simplex_vertex_range(f_simplex)) {
- simplex.insert(simplex.begin(), vertex);
- }
- filtrations.push_back(std::make_pair(simplex, Base::filtration(f_simplex)));
+ Simplex_and_filtration get_simplex_and_filtration(Simplex_handle f_simplex) {
+ Simplex simplex;
+ for (auto vertex : Base::simplex_vertex_range(f_simplex)) {
+ simplex.insert(simplex.begin(), vertex);
}
- return filtrations;
- }
-
- Filtered_simplices get_skeleton(int dimension) {
- Filtered_simplices skeletons;
- for (auto f_simplex : Base::skeleton_simplex_range(dimension)) {
- Simplex simplex;
- for (auto vertex : Base::simplex_vertex_range(f_simplex)) {
- simplex.insert(simplex.begin(), vertex);
- }
- skeletons.push_back(std::make_pair(simplex, Base::filtration(f_simplex)));
- }
- return skeletons;
+ return std::make_pair(std::move(simplex), Base::filtration(f_simplex));
}
Filtered_simplices get_star(const Simplex& simplex) {
@@ -131,9 +122,71 @@ 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
+ Complex_simplex_iterator get_simplices_iterator_begin() {
+ // this specific case works because the range is just a pair of iterators - won't work if range was a vector
+ return Base::complex_simplex_range().begin();
+ }
+
+ Complex_simplex_iterator get_simplices_iterator_end() {
+ // this specific case works because the range is just a pair of iterators - won't work if range was a vector
+ return Base::complex_simplex_range().end();
+ }
+
+ typename std::vector<Simplex_handle>::const_iterator get_filtration_iterator_begin() {
+ // Base::initialize_filtration(); already performed in filtration_simplex_range
+ // this specific case works because the range is just a pair of iterators - won't work if range was a vector
+ return Base::filtration_simplex_range().begin();
+ }
+
+ typename std::vector<Simplex_handle>::const_iterator get_filtration_iterator_end() {
+ // this specific case works because the range is just a pair of iterators - won't work if range was a vector
+ return Base::filtration_simplex_range().end();
+ }
+
+ Skeleton_simplex_iterator get_skeleton_iterator_begin(int dimension) {
+ // 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).begin();
+ }
+
+ Skeleton_simplex_iterator get_skeleton_iterator_end(int dimension) {
+ // 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();
}
};
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:
diff --git a/src/python/include/pybind11_diagram_utils.h b/src/python/include/pybind11_diagram_utils.h
new file mode 100644
index 00000000..d9627258
--- /dev/null
+++ b/src/python/include/pybind11_diagram_utils.h
@@ -0,0 +1,39 @@
+/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT.
+ * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details.
+ * Author(s): Marc Glisse
+ *
+ * Copyright (C) 2020 Inria
+ *
+ * Modification(s):
+ * - YYYY/MM Author: Description of the modification
+ */
+
+#include <pybind11/pybind11.h>
+#include <pybind11/numpy.h>
+
+#include <boost/range/counting_range.hpp>
+#include <boost/range/adaptor/transformed.hpp>
+
+namespace py = pybind11;
+typedef py::array_t<double> Dgm;
+
+// Get m[i,0] and m[i,1] as a pair
+static auto pairify(void* p, ssize_t h, ssize_t w) {
+ return [=](ssize_t i){
+ char* birth = (char*)p + i * h;
+ char* death = birth + w;
+ return std::make_pair(*(double*)birth, *(double*)death);
+ };
+}
+
+inline auto numpy_to_range_of_pairs(py::array_t<double> dgm) {
+ py::buffer_info buf = dgm.request();
+ // shape (n,2) or (0) for empty
+ if((buf.ndim!=2 || buf.shape[1]!=2) && (buf.ndim!=1 || buf.shape[0]!=0))
+ throw std::runtime_error("Diagram must be an array of size n x 2");
+ // In the case of shape (0), avoid reading non-existing strides[1] even if we won't use it.
+ ssize_t stride1 = buf.ndim == 2 ? buf.strides[1] : 0;
+ auto cnt = boost::counting_range<ssize_t>(0, buf.shape[0]);
+ return boost::adaptors::transform(cnt, pairify(buf.ptr, buf.strides[0], stride1));
+ // Be careful that the returned range cannot contain references to dead temporaries.
+}
generated by cgit v1.2.3 (git 2.39.1) at 2024-06-03 02:01:05 +0000