diff --git a/src/Persistent_cohomology/test/betti_numbers_unit_test.cpp b/src/Persistent_cohomology/test/betti_numbers_unit_test.cpp index eccd17b21f..96b880cc94 100644 --- a/src/Persistent_cohomology/test/betti_numbers_unit_test.cpp +++ b/src/Persistent_cohomology/test/betti_numbers_unit_test.cpp @@ -19,6 +19,8 @@ struct MiniSTOptions : Gudhi::Simplex_tree_options_full_featured { static const bool store_key = true; // I have few vertices typedef short Vertex_handle; + + static const bool is_multi_parameter = false; }; using Mini_simplex_tree = Gudhi::Simplex_tree; diff --git a/src/Persistent_cohomology/test/persistent_cohomology_unit_test.cpp b/src/Persistent_cohomology/test/persistent_cohomology_unit_test.cpp index 5a68ffb600..71c6717087 100644 --- a/src/Persistent_cohomology/test/persistent_cohomology_unit_test.cpp +++ b/src/Persistent_cohomology/test/persistent_cohomology_unit_test.cpp @@ -177,6 +177,7 @@ struct MiniSTOptions { static const bool contiguous_vertices = false; static const bool link_nodes_by_label = false; static const bool stable_simplex_handles = false; + static const bool is_multi_parameter = false; }; using Mini_simplex_tree = Gudhi::Simplex_tree; diff --git a/src/Simplex_tree/benchmark/simplex_tree_cofaces_benchmark.cpp b/src/Simplex_tree/benchmark/simplex_tree_cofaces_benchmark.cpp index f6e001b631..e4353977b6 100644 --- a/src/Simplex_tree/benchmark/simplex_tree_cofaces_benchmark.cpp +++ b/src/Simplex_tree/benchmark/simplex_tree_cofaces_benchmark.cpp @@ -91,6 +91,7 @@ struct Simplex_tree_options_stable_simplex_handles { static const bool contiguous_vertices = false; static const bool link_nodes_by_label = true; static const bool stable_simplex_handles = true; + static const bool is_multi_parameter = false; }; int main(int argc, char *argv[]) { diff --git a/src/Simplex_tree/concept/SimplexTreeOptions.h b/src/Simplex_tree/concept/SimplexTreeOptions.h index 0c95d3011c..45e5f03bca 100644 --- a/src/Simplex_tree/concept/SimplexTreeOptions.h +++ b/src/Simplex_tree/concept/SimplexTreeOptions.h @@ -31,5 +31,7 @@ struct SimplexTreeOptions { static const bool link_nodes_by_label; /// If true, Simplex_handle will not be invalidated after insertions or removals. static const bool stable_simplex_handles; + /// If true, assumes that Filtration_value is vector-like instead of float-like. This also assumes that Filtration_values is a class, which has a push_to method to push a filtration value $x$ onto $this>=0$. + static const bool is_multi_parameter; }; diff --git a/src/Simplex_tree/include/gudhi/Simplex_tree.h b/src/Simplex_tree/include/gudhi/Simplex_tree.h index ae4545876c..969148f0d6 100644 --- a/src/Simplex_tree/include/gudhi/Simplex_tree.h +++ b/src/Simplex_tree/include/gudhi/Simplex_tree.h @@ -154,16 +154,18 @@ class Simplex_tree { Key_simplex_base; struct Filtration_simplex_base_real { - Filtration_simplex_base_real() : filt_(0) {} - void assign_filtration(Filtration_value f) { filt_ = f; } - Filtration_value filtration() const { return filt_; } + Filtration_simplex_base_real() : filt_{} {} + void assign_filtration(const Filtration_value& f) { filt_ = f; } + const Filtration_value& filtration() const { return filt_; } + Filtration_value& filtration() { return filt_; } private: Filtration_value filt_; }; struct Filtration_simplex_base_dummy { Filtration_simplex_base_dummy() {} void assign_filtration(Filtration_value GUDHI_CHECK_code(f)) { GUDHI_CHECK(f == 0, "filtration value specified for a complex that does not store them"); } - Filtration_value filtration() const { return 0; } + const Filtration_value& filtration() const { return null_value; } + static constexpr Filtration_value null_value={}; }; typedef typename std::conditional::type Filtration_simplex_base; @@ -586,7 +588,7 @@ class Simplex_tree { * * Same as `filtration()`, but does not handle `null_simplex()`. */ - static Filtration_value filtration_(Simplex_handle sh) { + static const Filtration_value& filtration_(Simplex_handle sh) { GUDHI_CHECK (sh != null_simplex(), "null simplex"); return sh->second.filtration(); } @@ -614,18 +616,25 @@ class Simplex_tree { * Called on the null_simplex, it returns infinity. * If SimplexTreeOptions::store_filtration is false, returns 0. */ - static Filtration_value filtration(Simplex_handle sh) { + static const Filtration_value& filtration(Simplex_handle sh){ if (sh != null_simplex()) { return sh->second.filtration(); } else { - return std::numeric_limits::infinity(); + return inf_; + } + } + static Filtration_value& filtration_mutable(Simplex_handle sh){ + if (sh != null_simplex()) { + return sh->second.filtration(); + } else { + return inf_; } } /** \brief Sets the filtration value of a simplex. * \exception std::invalid_argument In debug mode, if sh is a null_simplex. */ - void assign_filtration(Simplex_handle sh, Filtration_value fv) { + void assign_filtration(Simplex_handle sh, const Filtration_value& fv) { GUDHI_CHECK(sh != null_simplex(), std::invalid_argument("Simplex_tree::assign_filtration - cannot assign filtration on null_simplex")); sh->second.assign_filtration(fv); @@ -848,7 +857,7 @@ class Simplex_tree { */ template > std::pair insert_simplex_raw(const RandomVertexHandleRange& simplex, - Filtration_value filtration) { + const Filtration_value& filtration) { Siblings * curr_sib = &root_; std::pair res_insert; auto vi = simplex.begin(); @@ -914,7 +923,7 @@ class Simplex_tree { * .end() return input iterators, with 'value_type' Vertex_handle. */ template> std::pair insert_simplex(const InputVertexRange & simplex, - Filtration_value filtration = 0) { + const Filtration_value& filtration = {}) { auto first = std::begin(simplex); auto last = std::end(simplex); @@ -943,7 +952,7 @@ class Simplex_tree { */ template> std::pair insert_simplex_and_subfaces(const InputVertexRange& Nsimplex, - Filtration_value filtration = 0) { + const Filtration_value& filtration = {}) { auto first = std::begin(Nsimplex); auto last = std::end(Nsimplex); @@ -972,7 +981,7 @@ class Simplex_tree { std::pair rec_insert_simplex_and_subfaces_sorted(Siblings* sib, ForwardVertexIterator first, ForwardVertexIterator last, - Filtration_value filt) { + const Filtration_value& filt) { // An alternative strategy would be: // - try to find the complete simplex, if found (and low filtration) exit // - insert all the vertices at once in sib @@ -989,8 +998,20 @@ class Simplex_tree { Simplex_handle simplex_one = insertion_result.first; bool one_is_new = insertion_result.second; if (!one_is_new) { - if (filtration(simplex_one) > filt) { - assign_filtration(simplex_one, filt); + if (!(filtration(simplex_one) <= filt)) { + if constexpr (SimplexTreeOptions::is_multi_parameter){ + // By default, does nothing and assumes that the user is smart. + if (filt < filtration(simplex_one)){ + // placeholder for comparable filtrations + } + else{ + // placeholder for incomparable filtrations + } + } + else{ // non-multiparameter + assign_filtration(simplex_one, filt); + } + } else { // FIXME: this interface makes no sense, and it doesn't seem to be tested. insertion_result.first = null_simplex(); @@ -1335,7 +1356,7 @@ class Simplex_tree { * The complex does not need to be empty before calling this function. However, if a vertex is * already present, its filtration value is not modified, unlike with other insertion functions. */ template - void insert_batch_vertices(VertexRange const& vertices, Filtration_value filt = 0) { + void insert_batch_vertices(VertexRange const& vertices, const Filtration_value& filt ={}) { auto verts = vertices | boost::adaptors::transformed([&](auto v){ return Dit_value_t(v, Node(&root_, filt)); }); root_.members_.insert(boost::begin(verts), boost::end(verts)); @@ -1689,7 +1710,7 @@ class Simplex_tree { static void intersection(std::vector >& intersection, Dictionary_it begin1, Dictionary_it end1, Dictionary_it begin2, Dictionary_it end2, - Filtration_value filtration_) { + const Filtration_value& filtration_) { if (begin1 == end1 || begin2 == end2) return; // ----->> while (true) { @@ -1901,12 +1922,22 @@ class Simplex_tree { if (dim == 0) return; // Find the maximum filtration value in the border Boundary_simplex_range&& boundary = boundary_simplex_range(sh); - Boundary_simplex_iterator max_border = std::max_element(std::begin(boundary), std::end(boundary), + Filtration_value max_filt_border_value; + if constexpr (SimplexTreeOptions::is_multi_parameter){ + // in that case, we assume that Filtration_value has a `push_to` member to handle this. + max_filt_border_value = Filtration_value(this->number_of_parameters_); + for (auto &face_sh : boundary){ + max_filt_border_value.push_to(filtration(face_sh)); // pushes the value of max_filt_border_value to reach simplex' filtration + } + } + else{ + Boundary_simplex_iterator max_border = std::max_element(std::begin(boundary), std::end(boundary), [](Simplex_handle sh1, Simplex_handle sh2) { return filtration(sh1) < filtration(sh2); }); - - Filtration_value max_filt_border_value = filtration(*max_border); + max_filt_border_value = filtration(*max_border); + } + // Replacing if(f=max)) would mean that if f is NaN, we replace it with the max of the children. // That seems more useful than keeping NaN. if (!(sh->second.filtration() >= max_filt_border_value)) { @@ -1941,7 +1972,7 @@ class Simplex_tree { * than it was before. However, `upper_bound_dimension()` will return the old value, which remains a valid upper * bound. If you care, you can call `dimension()` to recompute the exact dimension. */ - bool prune_above_filtration(Filtration_value filtration) { + bool prune_above_filtration(const Filtration_value& filtration) { if (std::numeric_limits::has_infinity && filtration == std::numeric_limits::infinity()) return false; // ---->> bool modified = rec_prune_above_filtration(root(), filtration); @@ -1951,7 +1982,7 @@ class Simplex_tree { } private: - bool rec_prune_above_filtration(Siblings* sib, Filtration_value filt) { + bool rec_prune_above_filtration(Siblings* sib, const Filtration_value& filt) { auto&& list = sib->members(); bool modified = false; bool emptied = false; @@ -2383,7 +2414,7 @@ class Simplex_tree { * @param[in] filt_value The new filtration value. * @param[in] min_dim The minimal dimension. Default value is 0. */ - void reset_filtration(Filtration_value filt_value, int min_dim = 0) { + void reset_filtration(const Filtration_value& filt_value, int min_dim = 0) { rec_reset_filtration(&root_, filt_value, min_dim); clear_filtration(); // Drop the cache. } @@ -2394,7 +2425,7 @@ class Simplex_tree { * @param[in] filt_value The new filtration value. * @param[in] min_depth The minimal depth. */ - void rec_reset_filtration(Siblings * sib, Filtration_value filt_value, int min_depth) { + void rec_reset_filtration(Siblings * sib, const Filtration_value& filt_value, int min_depth) { for (auto sh = sib->members().begin(); sh != sib->members().end(); ++sh) { if (min_depth <= 0) { sh->second.assign_filtration(filt_value); @@ -2560,6 +2591,18 @@ class Simplex_tree { /** \brief Upper bound on the dimension of the simplicial complex.*/ int dimension_; bool dimension_to_be_lowered_ = false; + +//MULTIPERS STUFF +public: + void set_number_of_parameters(int num){ + number_of_parameters_ = num; + } + int get_number_of_parameters() const{ + return number_of_parameters_; + } + inline static Filtration_value inf_ = std::numeric_limits::infinity(); +private: + int number_of_parameters_; }; // Print a Simplex_tree in os. @@ -2611,6 +2654,7 @@ struct Simplex_tree_options_full_featured { static const bool contiguous_vertices = false; static const bool link_nodes_by_label = false; static const bool stable_simplex_handles = false; + static const bool is_multi_parameter = false; }; /** Model of SimplexTreeOptions, faster than `Simplex_tree_options_full_featured` but note the unsafe @@ -2628,6 +2672,7 @@ struct Simplex_tree_options_fast_persistence { static const bool contiguous_vertices = true; static const bool link_nodes_by_label = false; static const bool stable_simplex_handles = false; + static const bool is_multi_parameter = false; }; /** Model of SimplexTreeOptions, faster cofaces than `Simplex_tree_options_full_featured`, note the @@ -2645,6 +2690,8 @@ struct Simplex_tree_options_fast_cofaces { static const bool contiguous_vertices = false; static const bool link_nodes_by_label = true; static const bool stable_simplex_handles = false; + static const bool is_multi_parameter = false; + }; /** @}*/ // end addtogroup simplex_tree diff --git a/src/Simplex_tree/include/gudhi/Simplex_tree/Simplex_tree_node_explicit_storage.h b/src/Simplex_tree/include/gudhi/Simplex_tree/Simplex_tree_node_explicit_storage.h index d7dbb541be..cb3205f816 100644 --- a/src/Simplex_tree/include/gudhi/Simplex_tree/Simplex_tree_node_explicit_storage.h +++ b/src/Simplex_tree/include/gudhi/Simplex_tree/Simplex_tree_node_explicit_storage.h @@ -41,7 +41,7 @@ struct GUDHI_EMPTY_BASE_CLASS_OPTIMIZATION Simplex_tree_node_explicit_storage : typedef typename SimplexTree::Simplex_key Simplex_key; Simplex_tree_node_explicit_storage(Siblings * sib = nullptr, - Filtration_value filtration = 0) + const Filtration_value& filtration = {}) : children_(sib) { this->assign_filtration(filtration); } diff --git a/src/Simplex_tree/test/simplex_tree_ctor_and_move_unit_test.cpp b/src/Simplex_tree/test/simplex_tree_ctor_and_move_unit_test.cpp index 986f811b54..e119008670 100644 --- a/src/Simplex_tree/test/simplex_tree_ctor_and_move_unit_test.cpp +++ b/src/Simplex_tree/test/simplex_tree_ctor_and_move_unit_test.cpp @@ -36,6 +36,7 @@ struct Simplex_tree_options_stable_simplex_handles { static const bool contiguous_vertices = false; static const bool link_nodes_by_label = true; static const bool stable_simplex_handles = true; + static const bool is_multi_parameter = false; }; typedef boost::mpl::list, diff --git a/src/Simplex_tree/test/simplex_tree_graph_expansion_unit_test.cpp b/src/Simplex_tree/test/simplex_tree_graph_expansion_unit_test.cpp index ed0f38d632..3a0eea485f 100644 --- a/src/Simplex_tree/test/simplex_tree_graph_expansion_unit_test.cpp +++ b/src/Simplex_tree/test/simplex_tree_graph_expansion_unit_test.cpp @@ -31,6 +31,7 @@ struct Simplex_tree_options_stable_simplex_handles { static const bool contiguous_vertices = false; static const bool link_nodes_by_label = false; static const bool stable_simplex_handles = true; + static const bool is_multi_parameter = false; }; typedef boost::mpl::list, diff --git a/src/Simplex_tree/test/simplex_tree_make_filtration_non_decreasing_unit_test.cpp b/src/Simplex_tree/test/simplex_tree_make_filtration_non_decreasing_unit_test.cpp index 05b7052612..3ac3b97401 100644 --- a/src/Simplex_tree/test/simplex_tree_make_filtration_non_decreasing_unit_test.cpp +++ b/src/Simplex_tree/test/simplex_tree_make_filtration_non_decreasing_unit_test.cpp @@ -33,6 +33,7 @@ struct Simplex_tree_options_stable_simplex_handles { static const bool contiguous_vertices = false; static const bool link_nodes_by_label = false; static const bool stable_simplex_handles = true; + static const bool is_multi_parameter = false; }; typedef boost::mpl::list, diff --git a/src/Simplex_tree/test/simplex_tree_unit_test.cpp b/src/Simplex_tree/test/simplex_tree_unit_test.cpp index fa1a444151..de1303e1e2 100644 --- a/src/Simplex_tree/test/simplex_tree_unit_test.cpp +++ b/src/Simplex_tree/test/simplex_tree_unit_test.cpp @@ -41,6 +41,7 @@ struct Simplex_tree_options_stable_simplex_handles { static const bool contiguous_vertices = false; static const bool link_nodes_by_label = true; static const bool stable_simplex_handles = true; + static const bool is_multi_parameter = false; }; typedef boost::mpl::list, diff --git a/src/python/CMakeLists.txt b/src/python/CMakeLists.txt index 721223cfd9..ac48f7eb39 100644 --- a/src/python/CMakeLists.txt +++ b/src/python/CMakeLists.txt @@ -55,6 +55,7 @@ if(PYTHONINTERP_FOUND) # Cython modules set(GUDHI_PYTHON_MODULES "${GUDHI_PYTHON_MODULES}'off_utils', ") set(GUDHI_PYTHON_MODULES "${GUDHI_PYTHON_MODULES}'simplex_tree', ") + set(GUDHI_PYTHON_MODULES "${GUDHI_PYTHON_MODULES}'simplex_tree_multi', ") set(GUDHI_PYTHON_MODULES "${GUDHI_PYTHON_MODULES}'edge_collapse', ") set(GUDHI_PYTHON_MODULES "${GUDHI_PYTHON_MODULES}'rips_complex', ") set(GUDHI_PYTHON_MODULES "${GUDHI_PYTHON_MODULES}'cubical_complex', ") @@ -165,6 +166,7 @@ if(PYTHONINTERP_FOUND) set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'off_utils', ") set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'simplex_tree', ") + set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'simplex_tree_multi', ") set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'rips_complex', ") set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'cubical_complex', ") set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'periodic_cubical_complex', ") diff --git a/src/python/gudhi/multiparameter_edge_collapse.py b/src/python/gudhi/multiparameter_edge_collapse.py new file mode 100644 index 0000000000..9337ce7321 --- /dev/null +++ b/src/python/gudhi/multiparameter_edge_collapse.py @@ -0,0 +1,29 @@ +from tqdm import tqdm +from filtration_domination import remove_strongly_filtration_dominated, remove_filtration_dominated + +def _collapse_edge_list(edges, num:int=0, full:bool=False, strong:bool=False, progress:bool=False): + """ + Given an edge list defining a 1 critical 2 parameter 1 dimensional simplicial complex, simplificates this filtered simplicial complex, using filtration-domination's edge collapser. + """ + n = len(edges) + if full: + num = 100 + with tqdm(range(num), total=num, desc="Removing edges", disable=not(progress)) as I: + for i in I: + if strong: + edges = remove_strongly_filtration_dominated(edges) # nogil ? + else: + edges = remove_filtration_dominated(edges) + # Prevents doing useless collapses + if len(edges) >= n: + if full and strong: + strong = False + n = len(edges) + # n = edges.size() # len(edges) + else : + break + else: + n = len(edges) + # n = edges.size() + return edges + diff --git a/src/python/gudhi/simplex_tree_multi.pxd b/src/python/gudhi/simplex_tree_multi.pxd new file mode 100644 index 0000000000..187633b1c9 --- /dev/null +++ b/src/python/gudhi/simplex_tree_multi.pxd @@ -0,0 +1,122 @@ +# 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): Vincent Rouvreau +# +# Copyright (C) 2016 Inria +# +# Modification(s): +# - 2022 David Loiseaux, Hannah Schreiber: adapt for multipersistence. +# - YYYY/MM Author: Description of the modification + +from cython cimport numeric +from libcpp.vector cimport vector +from libcpp.utility cimport pair +from libcpp cimport bool +from libcpp.string cimport string + +__author__ = "Vincent Rouvreau" +__copyright__ = "Copyright (C) 2016 Inria" +__license__ = "MIT" + +ctypedef int dimension_type +ctypedef float value_type +ctypedef Finitely_critical_multi_filtration filtration_type +ctypedef vector[value_type] python_filtration_type ## TODO: move constructor for C++ filtration type ? +ctypedef vector[int] simplex_type +ctypedef vector[simplex_type] simplex_list +ctypedef vector[pair[pair[int,int], pair[double, double]]] edge_list +ctypedef vector[int] euler_char_list +ctypedef pair[simplex_type, value_type*] simplex_filtration_type + + + + +cdef extern from "multi_filtrations/finitely_critical_filtrations.h" namespace "Gudhi::multiparameter::multi_filtrations": + cdef cppclass Finitely_critical_multi_filtration "Gudhi::multiparameter::multi_filtrations::Finitely_critical_multi_filtration": + Finitely_critical_multi_filtration() except + nogil + Finitely_critical_multi_filtration(vector[value_type]) except + + Finitely_critical_multi_filtration& operator=(const Finitely_critical_multi_filtration&) + filtration_type& get_vector() nogil + int size() nogil + void push_back(const value_type&) nogil + void clear() nogil + value_type& operator[](int) + + + +cdef extern from "Simplex_tree_interface_multi.h" namespace "Gudhi::multiparameter": + cdef cppclass Simplex_tree_options_multidimensional_filtration: + pass + + cdef cppclass Simplex_tree_multi_simplex_handle "Gudhi::multiparameter::Simplex_tree_interface_multi::Simplex_handle": + pass + + cdef cppclass Simplex_tree_multi_simplices_iterator "Gudhi::multiparameter::Simplex_tree_interface_multi::Complex_simplex_iterator": + Simplex_tree_multi_simplices_iterator() nogil + Simplex_tree_multi_simplex_handle& operator*() nogil + Simplex_tree_multi_simplices_iterator operator++() nogil + bint operator!=(Simplex_tree_multi_simplices_iterator) nogil + + cdef cppclass Simplex_tree_multi_skeleton_iterator "Gudhi::multiparameter::Simplex_tree_interface_multi::Skeleton_simplex_iterator": + Simplex_tree_multi_skeleton_iterator() nogil + Simplex_tree_multi_simplex_handle& operator*() nogil + Simplex_tree_multi_skeleton_iterator operator++() nogil + bint operator!=(Simplex_tree_multi_skeleton_iterator) nogil + + cdef cppclass Simplex_tree_multi_boundary_iterator "Gudhi::multiparameter::Simplex_tree_interface_multi::Boundary_simplex_iterator": + Simplex_tree_multi_boundary_iterator() nogil + Simplex_tree_multi_simplex_handle& operator*() nogil + Simplex_tree_multi_boundary_iterator operator++() nogil + bint operator!=(Simplex_tree_multi_boundary_iterator) nogil + + + cdef cppclass Simplex_tree_multi_interface "Gudhi::multiparameter::Simplex_tree_interface_multi": + Simplex_tree_multi_interface() nogil + Simplex_tree_multi_interface(Simplex_tree_multi_interface&) nogil + value_type* simplex_filtration(const vector[int]& simplex) nogil + void assign_simplex_filtration(vector[int]& simplex, const filtration_type& filtration) nogil + void initialize_filtration() nogil + int num_vertices() nogil + int num_simplices() nogil + void set_dimension(int dimension) nogil + dimension_type dimension() nogil + dimension_type upper_bound_dimension() nogil + bool find_simplex(vector[int]& simplex) nogil + bool insert(vector[int]& simplex, filtration_type& filtration) nogil + vector[simplex_filtration_type] get_star(const vector[int]& simplex) nogil + vector[simplex_filtration_type] get_cofaces(const vector[int]& simplex, int dimension) nogil + void expansion(int max_dim) except + nogil + void remove_maximal_simplex(simplex_type simplex) nogil + # bool prune_above_filtration(filtration_type filtration) nogil + bool prune_above_dimension(int dimension) nogil + bool make_filtration_non_decreasing() except + nogil + # void compute_extended_filtration() nogil + # Simplex_tree_multi_interface* collapse_edges(int nb_collapse_iteration) except + nogil + void reset_filtration(const filtration_type& filtration, int dimension) nogil + bint operator==(Simplex_tree_multi_interface) nogil + # Iterators over Simplex tree + simplex_filtration_type get_simplex_and_filtration(Simplex_tree_multi_simplex_handle f_simplex) nogil + Simplex_tree_multi_simplices_iterator get_simplices_iterator_begin() nogil + Simplex_tree_multi_simplices_iterator get_simplices_iterator_end() nogil + vector[Simplex_tree_multi_simplex_handle].const_iterator get_filtration_iterator_begin() nogil + vector[Simplex_tree_multi_simplex_handle].const_iterator get_filtration_iterator_end() nogil + Simplex_tree_multi_skeleton_iterator get_skeleton_iterator_begin(int dimension) nogil + Simplex_tree_multi_skeleton_iterator get_skeleton_iterator_end(int dimension) nogil + pair[Simplex_tree_multi_boundary_iterator, Simplex_tree_multi_boundary_iterator] get_boundary_iterators(vector[int] simplex) except + nogil + # Expansion with blockers + ctypedef bool (*blocker_func_t)(vector[int], void *user_data) + void expansion_with_blockers_callback(int dimension, blocker_func_t user_func, void *user_data) + + ## MULTIPERS STUFF + void set_keys_to_enumerate() nogil const + int get_key(const simplex_type) nogil + void set_key(simplex_type, int) nogil + void fill_lowerstar(const vector[value_type]&, int) nogil + simplex_list get_simplices_of_dimension(int) nogil + edge_list get_edge_list() nogil + # euler_char_list euler_char(const vector[filtration_type]&) nogil + void resize_all_filtrations(int) nogil + void set_number_of_parameters(int) nogil + int get_number_of_parameters() nogil + + diff --git a/src/python/gudhi/simplex_tree_multi.pyi b/src/python/gudhi/simplex_tree_multi.pyi new file mode 100644 index 0000000000..87b02aabd6 --- /dev/null +++ b/src/python/gudhi/simplex_tree_multi.pyi @@ -0,0 +1,715 @@ +import numpy as np +from gudhi.simplex_tree import SimplexTree ## Small hack for typing +from typing import Iterable +from tqdm import tqdm + + +# SimplexTree python interface +class SimplexTreeMulti: + """The simplex tree is an efficient and flexible data structure for + representing general (filtered) simplicial complexes. The data structure + is described in Jean-Daniel Boissonnat and Clément Maria. The Simplex + Tree: An Efficient Data Structure for General Simplicial Complexes. + Algorithmica, pages 1–22, 2014. + + This class is a multi-filtered, with keys, and non contiguous vertices version + of the simplex tree. + """ + + def __init__(self, other = None, num_parameters:int=2,default_values=[]): + """SimplexTreeMulti constructor. + + :param other: If `other` is `None` (default value), an empty `SimplexTreeMulti` is created. + If `other` is a `SimplexTree`, the `SimplexTreeMulti` is constructed from a deep copy of `other`. + If `other` is a `SimplexTreeMulti`, the `SimplexTreeMulti` is constructed from a deep copy of `other`. + :type other: SimplexTree or SimplexTreeMulti (Optional) + :param num_parameters: The number of parameter of the multi-parameter filtration. + :type num_parameters: int + :returns: An empty or a copy simplex tree. + :rtype: SimplexTreeMulti + + :raises TypeError: In case `other` is neither `None`, nor a `SimplexTree`, nor a `SimplexTreeMulti`. + """ + + + + def __is_defined(self): + """Returns true if SimplexTree pointer is not NULL. + """ + pass + + # def __is_persistence_defined(self): + # """Returns true if Persistence pointer is not NULL. + # """ + # return self.pcohptr != NULL + + def copy(self)->SimplexTreeMulti: + """ + :returns: A simplex tree that is a deep copy of itself. + :rtype: SimplexTreeMulti + + :note: The persistence information is not copied. If you need it in the clone, you have to call + :func:`compute_persistence` on it even if you had already computed it in the original. + """ + ... + + def __deepcopy__(self): + ... + + def filtration(self, simplex:list|np.ndarray)->np.ndarray: + """This function returns the filtration value for a given N-simplex in + this simplicial complex, or +infinity if it is not in the complex. + + :param simplex: The N-simplex, represented by a list of vertex. + :type simplex: list of int + :returns: The simplicial complex multi-critical filtration value. + :rtype: numpy array of shape (-1, num_parameters) + """ + ... + + def assign_filtration(self, simplex:list|np.ndarray, filtration:list|np.ndarray)->None: + """This function assigns a new multi-critical filtration value to a + given N-simplex. + + :param simplex: The N-simplex, represented by a list of vertex. + :type simplex: list of int + :param filtration: The new filtration(s) value(s), concatenated. + :type filtration: list[float] or np.ndarray[float, ndim=1] + + .. note:: + Beware that after this operation, the structure may not be a valid + filtration anymore, a simplex could have a lower filtration value + than one of its faces. Callers are responsible for fixing this + (with more :meth:`assign_filtration` or + :meth:`make_filtration_non_decreasing` for instance) before calling + any function that relies on the filtration property, like + :meth:`persistence`. + """ + ... + + def __getitem__(self, simplex): + ... + + + @property + def num_vertices(self)->int: + """This function returns the number of vertices of the simplicial + complex. + + :returns: The simplicial complex number of vertices. + :rtype: int + """ + ... + + @property + def num_simplices(self)->int: + """This function returns the number of simplices of the simplicial + complex. + + :returns: the simplicial complex number of simplices. + :rtype: int + """ + ... + + @property + def dimension(self)->int: + """This function returns the dimension of the simplicial complex. + + :returns: the simplicial complex dimension. + :rtype: int + + .. note:: + + This function is not constant time because it can recompute + dimension if required (can be triggered by + :func:`remove_maximal_simplex` + or + :func:`prune_above_filtration` + methods). + """ + ... + def upper_bound_dimension(self)->int: + """This function returns a valid dimension upper bound of the + simplicial complex. + + :returns: an upper bound on the dimension of the simplicial complex. + :rtype: int + """ + return self.get_ptr().upper_bound_dimension() + + def set_dimension(self, dimension)->None: + """This function sets the dimension of the simplicial complex. + + :param dimension: The new dimension value. + :type dimension: int + + .. note:: + + This function must be used with caution because it disables + dimension recomputation when required + (this recomputation can be triggered by + :func:`remove_maximal_simplex` + or + :func:`prune_above_filtration` + ). + """ + ... + + # def find(self, simplex)->bool: + # """This function returns if the N-simplex was found in the simplicial + # complex or not. + + # :param simplex: The N-simplex to find, represented by a list of vertex. + # :type simplex: list of int + # :returns: true if the simplex was found, false otherwise. + # :rtype: bool + # """ + # return self.get_ptr().find_simplex(simplex) + def __contains__(self, simplex)->bool: + """This function returns if the N-simplex was found in the simplicial + complex or not. + + :param simplex: The N-simplex to find, represented by a list of vertex. + :type simplex: list of int + :returns: true if the simplex was found, false otherwise. + :rtype: bool + """ + ... + + def insert(self, simplex, filtration:list|np.ndarray|None=None)->bool: + """This function inserts the given N-simplex and its subfaces with the + given filtration value (default value is '0.0'). If some of those + simplices are already present with a higher filtration value, their + filtration value is lowered. + + :param simplex: The N-simplex to insert, represented by a list of + vertex. + :type simplex: list of int + :param filtration: The filtration value of the simplex. + :type filtration: float + :returns: true if the simplex was not yet in the complex, false + otherwise (whatever its original filtration value). + :rtype: bool + """ + ... + + @cython.boundscheck(False) + @cython.wraparound(False) + def insert_batch(self, some_int[:,:] vertex_array, some_float[:,:] filtrations)->SimplexTreeMulti: + """Inserts k-simplices given by a sparse array in a format similar + to `torch.sparse `_. + The n-th simplex has vertices `vertex_array[0,n]`, ..., + `vertex_array[k,n]` and filtration value `filtrations[n,num_parameters]`. + /!\ Only compatible with 1-critical filtrations. If a simplex is repeated, + only one filtration value will be taken into account. + + :param vertex_array: the k-simplices to insert. + :type vertex_array: numpy.array of shape (k+1,n) + :param filtrations: the filtration values. + :type filtrations: numpy.array of shape (n,num_parameters) + """ + # TODO : multi-critical + # cdef vector[int] vertices = np.unique(vertex_array) + ... + + + @cython.boundscheck(False) + @cython.wraparound(False) + def assign_batch_filtration(self, some_int[:,:] vertex_array, some_float[:,:] filtrations, bool propagate=True)->SimplexTreeMulti: + """Assign k-simplices given by a sparse array in a format similar + to `torch.sparse `_. + The n-th simplex has vertices `vertex_array[0,n]`, ..., + `vertex_array[k,n]` and filtration value `filtrations[n,num_parameters]`. + /!\ Only compatible with 1-critical filtrations. If a simplex is repeated, + only one filtration value will be taken into account. + + :param vertex_array: the k-simplices to assign. + :type vertex_array: numpy.array of shape (k+1,n) + :param filtrations: the filtration values. + :type filtrations: numpy.array of shape (n,num_parameters) + """ + ... + + + + def get_simplices(self): + """This function returns a generator with simplices and their given + filtration values. + + :returns: The simplices. + :rtype: generator with tuples(simplex, filtration) + """ + ... + + def get_filtration(self): + """This function returns a generator with simplices and their given + filtration values sorted by increasing filtration values. + + :returns: The simplices sorted by increasing filtration values. + :rtype: generator with tuples(simplex, filtration) + """ + ... + + def get_skeleton(self, dimension): + """This function returns a generator with the (simplices of the) skeleton of a maximum given dimension. + + :param dimension: The skeleton dimension value. + :type dimension: int + :returns: The (simplices of the) skeleton of a maximum dimension. + :rtype: generator with tuples(simplex, filtration) + """ + ... + + def get_star(self, simplex): + """This function returns the star of a given N-simplex. + + :param simplex: The N-simplex, represented by a list of vertex. + :type simplex: list of int + :returns: The (simplices of the) star of a simplex. + :rtype: list of tuples(simplex, filtration) + """ + ... + + def get_cofaces(self, simplex, codimension): + """This function returns the cofaces of a given N-simplex with a + given codimension. + + :param simplex: The N-simplex, represented by a list of vertex. + :type simplex: list of int + :param codimension: The codimension. If codimension = 0, all cofaces + are returned (equivalent of get_star function) + :type codimension: int + :returns: The (simplices of the) cofaces of a simplex + :rtype: list of tuples(simplex, filtration) + """ + ... + + def get_boundaries(self, simplex): + """This function returns a generator with the boundaries of a given N-simplex. + If you do not need the filtration values, the boundary can also be obtained as + :code:`itertools.combinations(simplex,len(simplex)-1)`. + + :param simplex: The N-simplex, represented by a list of vertex. + :type simplex: list of int. + :returns: The (simplices of the) boundary of a simplex + :rtype: generator with tuples(simplex, filtration) + """ + ... + + def remove_maximal_simplex(self, simplex): + """This function removes a given maximal N-simplex from the simplicial + complex. + + :param simplex: The N-simplex, represented by a list of vertex. + :type simplex: list of int + + .. note:: + + The dimension of the simplicial complex may be lower after calling + remove_maximal_simplex than it was before. However, + :func:`upper_bound_dimension` + method will return the old value, which + remains a valid upper bound. If you care, you can call + :func:`dimension` + to recompute the exact dimension. + """ + ... + + # def prune_above_filtration(self, filtration)->bool: + # """Prune above filtration value given as parameter. + + # :param filtration: Maximum threshold value. + # :type filtration: float + # :returns: The filtration modification information. + # :rtype: bool + + + # .. note:: + + # Note that the dimension of the simplicial complex may be lower + # after calling + # :func:`prune_above_filtration` + # than it was before. However, + # :func:`upper_bound_dimension` + # will return the old value, which remains a + # valid upper bound. If you care, you can call + # :func:`dimension` + # method to recompute the exact dimension. + # """ + # return self.get_ptr().prune_above_filtration(filtration) + + def expansion(self, int max_dim)->SimplexTreeMulti: + """Expands the simplex tree containing only its one skeleton + until dimension max_dim. + + The expanded simplicial complex until dimension :math:`d` + attached to a graph :math:`G` is the maximal simplicial complex of + dimension at most :math:`d` admitting the graph :math:`G` as + :math:`1`-skeleton. + The filtration value assigned to a simplex is the maximal filtration + value of one of its edges. + + The simplex tree must contain no simplex of dimension bigger than + 1 when calling the method. + + :param max_dim: The maximal dimension. + :type max_dim: int + """ + ... + + def make_filtration_non_decreasing(self)->bool: + """This function ensures that each simplex has a higher filtration + value than its faces by increasing the filtration values. + + :returns: True if any filtration value was modified, + False if the filtration was already non-decreasing. + :rtype: bool + """ + ... + + def reset_filtration(self, filtration, min_dim = 0): + """This function resets the filtration value of all the simplices of dimension at least min_dim. Resets all the + simplex tree when `min_dim = 0`. + `reset_filtration` may break the filtration property with `min_dim > 0`, and it is the user's responsibility to + make it a valid filtration (using a large enough `filt_value`, or calling `make_filtration_non_decreasing` + afterwards for instance). + + :param filtration: New threshold value. + :type filtration: float. + :param min_dim: The minimal dimension. Default value is 0. + :type min_dim: int. + """ + ... + + + + # def extend_filtration(self): + # """ Extend filtration for computing extended persistence. This function only uses the filtration values at the + # 0-dimensional simplices, and computes the extended persistence diagram induced by the lower-star filtration + # computed with these values. + # + # .. note:: + # + # Note that after calling this function, the filtration values are actually modified within the simplex tree. + # The function :func:`extended_persistence` retrieves the original values. + # + # .. note:: + # + # Note that this code creates an extra vertex internally, so you should make sure that the simplex tree does + # not contain a vertex with the largest possible value (i.e., 4294967295). + # + # This `notebook `_ + # explains how to compute an extension of persistence called extended persistence. + # """ + # self.get_ptr().compute_extended_filtration() + + # def extended_persistence(self, homology_coeff_field=11, min_persistence=0): + # """This function retrieves good values for extended persistence, and separate the diagrams into the Ordinary, + # Relative, Extended+ and Extended- subdiagrams. + # + # :param homology_coeff_field: The homology coefficient field. Must be a prime number. Default value is 11. Max is 46337. + # :type homology_coeff_field: int + # :param min_persistence: The minimum persistence value (i.e., the absolute value of the difference between the + # persistence diagram point coordinates) to take into account (strictly greater than min_persistence). + # Default value is 0.0. Sets min_persistence to -1.0 to see all values. + # :type min_persistence: float + # :returns: A list of four persistence diagrams in the format described in :func:`persistence`. The first one is + # Ordinary, the second one is Relative, the third one is Extended+ and the fourth one is Extended-. + # See https://link.springer.com/article/10.1007/s10208-008-9027-z and/or section 2.2 in + # https://link.springer.com/article/10.1007/s10208-017-9370-z for a description of these subtypes. + # + # .. note:: + # + # This function should be called only if :func:`extend_filtration` has been called first! + # + # .. note:: + # + # The coordinates of the persistence diagram points might be a little different than the + # original filtration values due to the internal transformation (scaling to [-2,-1]) that is + # performed on these values during the computation of extended persistence. + # + # This `notebook `_ + # explains how to compute an extension of persistence called extended persistence. + # """ + # cdef vector[pair[int, pair[value_type, value_type]]] persistence_result + # if self.pcohptr != NULL: + # del self.pcohptr + # self.pcohptr = new Simplex_tree_persistence_interface(self.get_ptr(), False) + # self.pcohptr.compute_persistence(homology_coeff_field, -1.) + # return self.pcohptr.compute_extended_persistence_subdiagrams(min_persistence) + + # TODO : cython3 + # def expansion_with_blocker(self, max_dim, blocker_func): + # """Expands the Simplex_tree containing only a graph. Simplices corresponding to cliques in the graph are added + # incrementally, faces before cofaces, unless the simplex has dimension larger than `max_dim` or `blocker_func` + # returns `True` for this simplex. + + # The function identifies a candidate simplex whose faces are all already in the complex, inserts it with a + # filtration value corresponding to the maximum of the filtration values of the faces, then calls `blocker_func` + # with this new simplex (represented as a list of int). If `blocker_func` returns `True`, the simplex is removed, + # otherwise it is kept. The algorithm then proceeds with the next candidate. + + # .. warning:: + # Several candidates of the same dimension may be inserted simultaneously before calling `blocker_func`, so + # if you examine the complex in `blocker_func`, you may hit a few simplices of the same dimension that have + # not been vetted by `blocker_func` yet, or have already been rejected but not yet removed. + + # :param max_dim: Expansion maximal dimension value. + # :type max_dim: int + # :param blocker_func: Blocker oracle. + # :type blocker_func: Callable[[List[int]], bool] + # """ + # self.get_ptr().expansion_with_blockers_callback(max_dim, callback, blocker_func) + + # def persistence(self, homology_coeff_field=11, min_persistence=0, persistence_dim_max = False): + # """This function computes and returns the persistence of the simplicial complex. + # + # :param homology_coeff_field: The homology coefficient field. Must be a + # prime number. Default value is 11. Max is 46337. + # :type homology_coeff_field: int + # :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 + # :param persistence_dim_max: If true, the persistent homology for the + # maximal dimension in the complex is computed. If false, it is + # ignored. Default is false. + # :type persistence_dim_max: bool + # :returns: The persistence of the simplicial complex. + # :rtype: list of pairs(dimension, pair(birth, death)) + # """ + # self.compute_persistence(homology_coeff_field, min_persistence, persistence_dim_max) + # return self.pcohptr.get_persistence() + + + +## This function is only meant for the edge collapse interface. + def get_edge_list(self): + ... + + def collapse_edges(self, max_dimension:int=None, num:int=1, progress:bool=False, strong:bool=True, full:bool=False, ignore_warning:bool=False)->SimplexTreeMulti: + """Edge collapse for 1-critical 2-parameter clique complex (see https://arxiv.org/abs/2211.05574). + It uses the code from the github repository https://github.com/aj-alonso/filtration_domination . + + Parameters + ---------- + max_dimension:int + Max simplicial dimension of the complex. Unless specified, keeps the same dimension. + num:int + The number of collapses to do. + strong:bool + Whether to use strong collapses or standard collapses (slower, but may remove more edges) + full:bool + Collapses the maximum number of edges if true, i.e., will do (at most) 100 strong collapses and (at most) 100 non-strong collapses afterward. + progress:bool + If true, shows the progress of the number of collapses. + + WARNING + ------- + - This will destroy all of the k-simplices, with k>=2. Be sure to use this with a clique complex, if you want to preserve the homology >= dimension 1. + - This is for 1 critical simplices, with 2 parameter persistence. + Returns + ------- + self:SimplexTreeMulti + A (smaller) simplex tree that has the same homology over this bifiltration. + + """ + # TODO : find a way to do multiple edge collapses without python conversions. + ... + + def _reconstruct_from_edge_list(self, edges, swap:bool=True, expand_dimension:int=None)->SimplexTreeMulti: + """ + Generates a 1-dimensional copy of self, with the edges given as input. Useful for edge collapses + + Input + ----- + - edges : Iterable[(int,int),(float,float)] ## This is the format of the rust library filtration-domination + - swap : bool + If true, will swap self and the collapsed simplextrees. + - expand_dim : int + expands back the simplextree to this dimension + Ouput + ----- + The reduced SimplexTreeMulti having only these edges. + """ + ... + + @property + def num_parameters(self)->int: + ... + def get_simplices_of_dimension(self, dim:int)->np.ndarray: + ... + def key(self, simplex:list|np.ndarray): + ... + def set_keys_to_enumerate(self)->None: + ... + def set_key(self,simplex:list|np.ndarray, key:int)->None: + ... + + + def to_scc(self, path="scc_dataset.txt", progress:bool=True, overwrite:bool=False, ignore_last_generators:bool=True, strip_comments:bool=False, reverse_block:bool=True, rivet_compatible=False)->None: + """ Create a file with the scc2020 standard, representing the n-filtration of the simplextree. + Link : https://bitbucket.org/mkerber/chain_complex_format/src/master/ + + Parameters + ---------- + path:str + path of the file. + ignore_last_generators:bool = True + If false, will include the filtration values of the last free persistence module. + progress:bool = True + Shows the progress bar. + overwrite:bool = False + If true, will overwrite the previous file if it already exists. + ignore_last_generators:bool=True + If true, does not write the final generators to the file. Rivet ignores them. + reverse_block:bool=True + Some obscure programs reverse the inside-block order. + rivet_compatible:bool=False + Returns a firep (old scc2020) format instead. Only Rivet uses this. + + Returns + ------- + Nothing + """ + ... + + def to_rivet(self, path="rivet_dataset.txt", degree:int|None = None, progress:bool=False, overwrite:bool=False, xbins:int|None=None, ybins:int|None=None)->None: + """ Create a file that can be imported by rivet, representing the filtration of the simplextree. + + Parameters + ---------- + path:str + path of the file. + degree:int + The homological degree to ask rivet to compute. + progress:bool = True + Shows the progress bar. + overwrite:bool = False + If true, will overwrite the previous file if it already exists. + Returns + ------- + Nothing + """ + ... + + + + def _get_filtration_values(self, vector[int] degrees, bool inf_to_nan:bool=False)->Iterable[np.ndarray]: + # cdef vector[int] c_degrees = degrees + ... + + @staticmethod + def _reduce_grid(filtrations_values,resolutions=None, strategy:str="exact", bool unique=True, some_float _q_factor=1., drop_quantiles=[0,0]): + ... + + def get_filtration_grid(self, resolution:Iterable[int]|None=None, degrees:Iterable[int]|None=None, drop_quantiles:float|tuple=0, grid_strategy:str="exact")->Iterable[np.ndarray]: + """ + Returns a grid over the n-filtration, from the simplextree. Usefull for grid_squeeze. TODO : multicritical + + Parameters + ---------- + resolution: list[int] + resolution of the grid, for each parameter + box=None : pair[list[float]] + Grid bounds. format : [low bound, high bound] + If None is given, will use the filtration bounds of the simplextree. + grid_strategy="regular" : string + Either "regular", "quantile", or "exact". + Returns + ------- + List of filtration values, for each parameter, defining the grid. + """ + ... + + + + def grid_squeeze(self, filtration_grid:np.ndarray|list|None=None, coordinate_values:bool=True, force=False, **filtration_grid_kwargs)->SimplexTreeMulti: + """ + Fit the filtration of the simplextree to a grid. + + :param filtration_grid: The grid on which to squeeze. An example of grid can be given by the `get_filtration_grid` method. + :type filtration_grid: list[list[float]] + :param coordinate_values: If true, the filtrations values of the simplices will be set to the coordinate of the filtration grid. + :type coordinate_values: bool + """ + ... + + @property + def _is_squeezed(self)->bool: + ... + + def filtration_bounds(self, degrees:Iterable[int]|None=None, q:float|tuple=0, split_dimension:bool=False)->np.ndarray: + """ + Returns the filtrations bounds of the finite filtration values. + """ + ... + + + + + def fill_lowerstar(self, F, parameter:int)->SimplexTreeMulti: + """ Fills the `dimension`th filtration by the lower-star filtration defined by F. + + Parameters + ---------- + F:1d array + The density over the vertices, that induces a lowerstar filtration. + parameter:int + Which filtration parameter to fill. /!\ python starts at 0. + + Returns + ------- + self:SimplexTreeMulti + """ + # for s, sf in self.get_simplices(): + # self.assign_filtration(s, [f if i != dimension else np.max(np.array(F)[s]) for i,f in enumerate(sf)]) + ... + + def project_on_line(self, parameter:int=0, basepoint:None|list|np.ndarray= None)->SimplexTree: + """Converts an multi simplextree to a gudhi simplextree. + Parameters + ---------- + parameter:int = 0 + The parameter to keep. WARNING will crash if the multi simplextree is not well filled. + basepoint:None + Instead of keeping a single parameter, will consider the filtration defined by the diagonal line crossing the basepoint. + WARNING + ------- + There are no safeguard yet, it WILL crash if asking for a parameter that is not filled. + Returns + ------- + A SimplexTree with chosen 1D filtration. + """ + ... + + def linear_projections(self, linear_forms:np.ndarray)->Iterable[SimplexTree]: + """ + Compute the 1-parameter projections, w.r.t. given the linear forms, of this simplextree. + + Input + ----- + - Array of shape (num_linear_forms, num_parameters) + + Output + ------ + - List of projected (gudhi) simplextrees. + """ + ... + + + def set_num_parameter(self, num:int): + """ + Sets the numbers of parameters. + WARNING : it will resize all the filtrations to this size. + """ + ... + + def __eq__(self, other:SimplexTreeMulti): + """Test for structural equality + :returns: True if the 2 simplex trees are equal, False otherwise. + :rtype: bool + """ + ... + \ No newline at end of file diff --git a/src/python/gudhi/simplex_tree_multi.pyx b/src/python/gudhi/simplex_tree_multi.pyx new file mode 100644 index 0000000000..d2beda7069 --- /dev/null +++ b/src/python/gudhi/simplex_tree_multi.pyx @@ -0,0 +1,1240 @@ +# 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): Vincent Rouvreau +# +# Copyright (C) 2016 Inria +# +# Modification(s): +# - 2023 David Loiseaux : Conversions with standard simplextree, scc2020 format, edge collapses, euler characteristic, grid filtrations. +# - 2022/11 Hannah Schreiber / David Loiseaux : adapt for multipersistence. +# - YYYY/MM Author: Description of the modification + + + +__author__ = "Vincent Rouvreau" +__copyright__ = "Copyright (C) 2016 Inria" +__license__ = "MIT" + +from libc.stdint cimport intptr_t, int32_t, int64_t +from cython.operator import dereference, preincrement +from libc.stdint cimport intptr_t +from libc.stdint cimport uintptr_t +from libcpp.map cimport map + + +ctypedef fused some_int: + int32_t + int64_t + int + +ctypedef fused some_float: + float + double + +ctypedef vector[pair[pair[int,int],pair[value_type,value_type]]] edge_list_type + +from typing import Any + +cimport numpy as cnp +import numpy as np +cnp.import_array() + +from gudhi.simplex_tree_multi cimport * +cimport cython +from gudhi.simplex_tree import SimplexTree ## Small hack for typing +from typing import Iterable +from tqdm import tqdm + + + +from warnings import warn + + +cdef extern from "Simplex_tree_multi.h" namespace "Gudhi::multiparameter": + void multify_from_ptr(const uintptr_t, const uintptr_t, const unsigned int, const vector[value_type]&) except + nogil + void flatten_from_ptr(const uintptr_t, const uintptr_t, const unsigned int) nogil + void linear_projection_from_ptr(const uintptr_t, const uintptr_t, const vector[value_type]&) nogil + void flatten_diag_from_ptr(const uintptr_t, const uintptr_t, const vector[value_type], int) nogil + void squeeze_filtration(uintptr_t, const vector[vector[value_type]]&, bool) except + nogil + vector[vector[vector[value_type]]] get_filtration_values(uintptr_t, const vector[int]&) except + nogil + + +# cdef bool callback(vector[int] simplex, void *blocker_func): + # return (blocker_func)(simplex) + +# SimplexTree python interface +cdef class SimplexTreeMulti: + """The simplex tree is an efficient and flexible data structure for + representing general (filtered) simplicial complexes. The data structure + is described in Jean-Daniel Boissonnat and Clément Maria. The Simplex + Tree: An Efficient Data Structure for General Simplicial Complexes. + Algorithmica, pages 1–22, 2014. + + This class is a multi-filtered, with keys, and non contiguous vertices version + of the simplex tree. + """ + # unfortunately 'cdef public Simplex_tree_multi_interface* thisptr' is not possible + # Use intptr_t instead to cast the pointer + cdef public intptr_t thisptr + cdef public vector[vector[value_type]] filtration_grid + + # Get the pointer casted as it should be + cdef Simplex_tree_multi_interface* get_ptr(self) nogil: + return (self.thisptr) + + # cdef Simplex_tree_persistence_interface * pcohptr + # Fake constructor that does nothing but documenting the constructor + def __init__(self, other = None, num_parameters:int=2,default_values=[], safe_conversion=False): + """SimplexTreeMulti constructor. + + :param other: If `other` is `None` (default value), an empty `SimplexTreeMulti` is created. + If `other` is a `SimplexTree`, the `SimplexTreeMulti` is constructed from a deep copy of `other`. + If `other` is a `SimplexTreeMulti`, the `SimplexTreeMulti` is constructed from a deep copy of `other`. + :type other: SimplexTree or SimplexTreeMulti (Optional) + :param num_parameters: The number of parameter of the multi-parameter filtration. + :type num_parameters: int + :returns: An empty or a copy simplex tree. + :rtype: SimplexTreeMulti + + :raises TypeError: In case `other` is neither `None`, nor a `SimplexTree`, nor a `SimplexTreeMulti`. + """ + + # The real cython constructor + def __cinit__(self, other = None, int num_parameters=2, + default_values=[-np.inf], # I'm not sure why `[]` does not work. Cython bug ? +bool safe_conversion=False, + ): #TODO doc + cdef vector[value_type] c_default_values=default_values + if other is not None: + if isinstance(other, SimplexTreeMulti): + self.thisptr = _get_copy_intptr(other) + num_parameters = other.num_parameters + self.filtration_grid = other.filtration_grid + elif isinstance(other, SimplexTree): # Constructs a SimplexTreeMulti from a SimplexTree + self.thisptr = (new Simplex_tree_multi_interface()) + if safe_conversion: + new_st_multi = _safe_simplextree_multify(other, num_parameters = num_parameters) + self.thisptr, new_st_multi.thisptr = new_st_multi.thisptr, self.thisptr + else: + multify_from_ptr(other.thisptr, self.thisptr, num_parameters, c_default_values) + else: + raise TypeError("`other` argument requires to be of type `SimplexTree`, `SimplexTreeMulti`, or `None`.") + else: + self.thisptr = (new Simplex_tree_multi_interface()) + self.get_ptr().set_number_of_parameters(num_parameters) + self.filtration_grid=[[]*num_parameters] + + def __dealloc__(self): + cdef Simplex_tree_multi_interface* ptr = self.get_ptr() + if ptr != NULL: + del ptr + # if self.pcohptr != NULL: + # del self.pcohptr + + def __is_defined(self): + """Returns true if SimplexTree pointer is not NULL. + """ + return self.get_ptr() != NULL + + # def __is_persistence_defined(self): + # """Returns true if Persistence pointer is not NULL. + # """ + # return self.pcohptr != NULL + + def copy(self)->SimplexTreeMulti: + """ + :returns: A simplex tree that is a deep copy of itself. + :rtype: SimplexTreeMulti + + :note: The persistence information is not copied. If you need it in the clone, you have to call + :func:`compute_persistence` on it even if you had already computed it in the original. + """ + stree = SimplexTreeMulti(self,num_parameters=self.num_parameters) + return stree + + def __deepcopy__(self): + return self.copy() + + def filtration(self, simplex:list|np.ndarray)->np.ndarray: + """This function returns the filtration value for a given N-simplex in + this simplicial complex, or +infinity if it is not in the complex. + + :param simplex: The N-simplex, represented by a list of vertex. + :type simplex: list of int + :returns: The simplicial complex multi-critical filtration value. + :rtype: numpy array of shape (-1, num_parameters) + """ + return self[simplex] + + def assign_filtration(self, simplex:list|np.ndarray, filtration:list|np.ndarray)->None: + """This function assigns a new multi-critical filtration value to a + given N-simplex. + + :param simplex: The N-simplex, represented by a list of vertex. + :type simplex: list of int + :param filtration: The new filtration(s) value(s), concatenated. + :type filtration: list[float] or np.ndarray[float, ndim=1] + + .. note:: + Beware that after this operation, the structure may not be a valid + filtration anymore, a simplex could have a lower filtration value + than one of its faces. Callers are responsible for fixing this + (with more :meth:`assign_filtration` or + :meth:`make_filtration_non_decreasing` for instance) before calling + any function that relies on the filtration property, like + :meth:`persistence`. + """ + assert len(filtration)>0 and len(filtration) % self.get_ptr().get_number_of_parameters() == 0 + self.get_ptr().assign_simplex_filtration(simplex, Finitely_critical_multi_filtration(filtration)) + + def __getitem__(self, simplex): + cdef vector[int] csimplex = simplex + cdef value_type[:] filtration_view = self.get_ptr().simplex_filtration(csimplex) + return np.asarray(filtration_view) + + + @property + def num_vertices(self)->int: + """This function returns the number of vertices of the simplicial + complex. + + :returns: The simplicial complex number of vertices. + :rtype: int + """ + return self.get_ptr().num_vertices() + + @property + def num_simplices(self)->int: + """This function returns the number of simplices of the simplicial + complex. + + :returns: the simplicial complex number of simplices. + :rtype: int + """ + return self.get_ptr().num_simplices() + + @property + def dimension(self)->int: + """This function returns the dimension of the simplicial complex. + + :returns: the simplicial complex dimension. + :rtype: int + + .. note:: + + This function is not constant time because it can recompute + dimension if required (can be triggered by + :func:`remove_maximal_simplex` + or + :func:`prune_above_filtration` + methods). + """ + return self.get_ptr().dimension() + def upper_bound_dimension(self)->int: + """This function returns a valid dimension upper bound of the + simplicial complex. + + :returns: an upper bound on the dimension of the simplicial complex. + :rtype: int + """ + return self.get_ptr().upper_bound_dimension() + + def set_dimension(self, dimension)->None: + """This function sets the dimension of the simplicial complex. + + :param dimension: The new dimension value. + :type dimension: int + + .. note:: + + This function must be used with caution because it disables + dimension recomputation when required + (this recomputation can be triggered by + :func:`remove_maximal_simplex` + or + :func:`prune_above_filtration` + ). + """ + self.get_ptr().set_dimension(dimension) + + # def find(self, simplex)->bool: + # """This function returns if the N-simplex was found in the simplicial + # complex or not. + + # :param simplex: The N-simplex to find, represented by a list of vertex. + # :type simplex: list of int + # :returns: true if the simplex was found, false otherwise. + # :rtype: bool + # """ + # return self.get_ptr().find_simplex(simplex) + def __contains__(self, simplex)->bool: + """This function returns if the N-simplex was found in the simplicial + complex or not. + + :param simplex: The N-simplex to find, represented by a list of vertex. + :type simplex: list of int + :returns: true if the simplex was found, false otherwise. + :rtype: bool + """ + return self.get_ptr().find_simplex(simplex) + + def insert(self, simplex, filtration:list|np.ndarray|None=None)->bool: + """This function inserts the given N-simplex and its subfaces with the + given filtration value (default value is '0.0'). If some of those + simplices are already present with a higher filtration value, their + filtration value is lowered. + + :param simplex: The N-simplex to insert, represented by a list of + vertex. + :type simplex: list of int + :param filtration: The filtration value of the simplex. + :type filtration: float + :returns: true if the simplex was not yet in the complex, false + otherwise (whatever its original filtration value). + :rtype: bool + """ + # TODO C++, to be compatible with insert_batch and multicritical filtrations + num_parameters = self.get_ptr().get_number_of_parameters() + assert filtration is None or len(filtration) % num_parameters == 0 + if filtration is None: + filtration = np.array([-np.inf]*num_parameters, dtype = float) + return self.get_ptr().insert(simplex, Finitely_critical_multi_filtration(filtration)) + + @cython.boundscheck(False) + @cython.wraparound(False) + def insert_batch(self, some_int[:,:] vertex_array, some_float[:,:] filtrations)->SimplexTreeMulti: + """Inserts k-simplices given by a sparse array in a format similar + to `torch.sparse `_. + The n-th simplex has vertices `vertex_array[0,n]`, ..., + `vertex_array[k,n]` and filtration value `filtrations[n,num_parameters]`. + /!\ Only compatible with 1-critical filtrations. If a simplex is repeated, + only one filtration value will be taken into account. + + :param vertex_array: the k-simplices to insert. + :type vertex_array: numpy.array of shape (k+1,n) + :param filtrations: the filtration values. + :type filtrations: numpy.array of shape (n,num_parameters) + """ + # TODO : multi-critical + # cdef vector[int] vertices = np.unique(vertex_array) + cdef Py_ssize_t k = vertex_array.shape[0] + cdef Py_ssize_t n = vertex_array.shape[1] + assert filtrations.shape[0] == n, 'inconsistent sizes for vertex_array and filtrations' + assert filtrations.shape[1] == self.num_parameters, "wrong number of parameters" + cdef Py_ssize_t i + cdef Py_ssize_t j + cdef vector[int] v + cdef Finitely_critical_multi_filtration w + cdef int n_parameters = self.num_parameters + with nogil: + for i in range(n): + for j in range(k): + v.push_back(vertex_array[j, i]) + for j in range(n_parameters): + w.push_back(filtrations[i,j]) + self.get_ptr().insert(v, w) + v.clear() + w.clear() + return self + + + @cython.boundscheck(False) + @cython.wraparound(False) + def assign_batch_filtration(self, some_int[:,:] vertex_array, some_float[:,:] filtrations, bool propagate=True)->SimplexTreeMulti: + """Assign k-simplices given by a sparse array in a format similar + to `torch.sparse `_. + The n-th simplex has vertices `vertex_array[0,n]`, ..., + `vertex_array[k,n]` and filtration value `filtrations[n,num_parameters]`. + /!\ Only compatible with 1-critical filtrations. If a simplex is repeated, + only one filtration value will be taken into account. + + :param vertex_array: the k-simplices to assign. + :type vertex_array: numpy.array of shape (k+1,n) + :param filtrations: the filtration values. + :type filtrations: numpy.array of shape (n,num_parameters) + """ + cdef Py_ssize_t k = vertex_array.shape[0] + cdef Py_ssize_t n = vertex_array.shape[1] + assert filtrations.shape[0] == n, 'inconsistent sizes for vertex_array and filtrations' + assert filtrations.shape[1] == self.num_parameters, "wrong number of parameters" + cdef Py_ssize_t i + cdef Py_ssize_t j + cdef vector[int] v + cdef Finitely_critical_multi_filtration w + cdef int n_parameters = self.num_parameters + with nogil: + for i in range(n): + for j in range(k): + v.push_back(vertex_array[j, i]) + for j in range(n_parameters): + w.push_back(filtrations[i,j]) + self.get_ptr().assign_simplex_filtration(v, w) + v.clear() + w.clear() + if propagate: self.make_filtration_non_decreasing() + return self + + + + def get_simplices(self): + """This function returns a generator with simplices and their given + filtration values. + + :returns: The simplices. + :rtype: generator with tuples(simplex, filtration) + """ + cdef Simplex_tree_multi_simplices_iterator it = self.get_ptr().get_simplices_iterator_begin() + cdef Simplex_tree_multi_simplices_iterator end = self.get_ptr().get_simplices_iterator_end() + cdef Simplex_tree_multi_simplex_handle sh = dereference(it) + # cdef pair[simplex_type,Finitely_critical_multi_filtration] out_ + # while it != end: + # out_ = self.get_ptr().get_simplex_and_filtration(dereference(it)) + # out = (out_.first,out_.second.get_vector()) + # yield out + # preincrement(it) + # cdef pair[simplex_type,filtration_type] out + cdef int num_parameters = self.get_ptr().get_number_of_parameters() + while it != end: + pair = self.get_ptr().get_simplex_and_filtration(dereference(it)) + + yield (np.asarray(pair.first, dtype=int),np.asarray( pair.second)) + # yield SimplexTreeMulti._pair_simplex_filtration_to_python(out) + preincrement(it) + + def get_filtration(self): + """This function returns a generator with simplices and their given + filtration values sorted by increasing filtration values. + + :returns: The simplices sorted by increasing filtration values. + :rtype: generator with tuples(simplex, filtration) + """ + cdef vector[Simplex_tree_multi_simplex_handle].const_iterator it = self.get_ptr().get_filtration_iterator_begin() + cdef vector[Simplex_tree_multi_simplex_handle].const_iterator end = self.get_ptr().get_filtration_iterator_end() + cdef int num_parameters = self.get_ptr().get_number_of_parameters() + while it != end: + # yield self.get_ptr().get_simplex_and_filtration(dereference(it)) + pair = self.get_ptr().get_simplex_and_filtration(dereference(it)) + yield (np.asarray(pair.first, dtype=int),np.asarray( pair.second)) + preincrement(it) + + def get_skeleton(self, dimension): + """This function returns a generator with the (simplices of the) skeleton of a maximum given dimension. + + :param dimension: The skeleton dimension value. + :type dimension: int + :returns: The (simplices of the) skeleton of a maximum dimension. + :rtype: generator with tuples(simplex, filtration) + """ + cdef Simplex_tree_multi_skeleton_iterator it = self.get_ptr().get_skeleton_iterator_begin(dimension) + cdef Simplex_tree_multi_skeleton_iterator end = self.get_ptr().get_skeleton_iterator_end(dimension) + cdef int num_parameters = self.get_ptr().get_number_of_parameters() + while it != end: + # yield self.get_ptr().get_simplex_and_filtration(dereference(it)) + pair = self.get_ptr().get_simplex_and_filtration(dereference(it)) + yield (np.asarray(pair.first, dtype=int),np.asarray( pair.second)) + preincrement(it) + + def get_star(self, simplex): + """This function returns the star of a given N-simplex. + + :param simplex: The N-simplex, represented by a list of vertex. + :type simplex: list of int + :returns: The (simplices of the) star of a simplex. + :rtype: list of tuples(simplex, filtration) + """ + cdef simplex_type csimplex = simplex + cdef int num_parameters = self.num_parameters + # for i in simplex: + # csimplex.push_back(i) + cdef vector[simplex_filtration_type] star \ + = self.get_ptr().get_star(csimplex) + ct = [] + + for filtered_simplex in star: + v = [] + for vertex in filtered_simplex.first: + v.append(vertex) + ct.append((v, np.asarray(filtered_simplex.second))) + return ct + + def get_cofaces(self, simplex, codimension): + """This function returns the cofaces of a given N-simplex with a + given codimension. + + :param simplex: The N-simplex, represented by a list of vertex. + :type simplex: list of int + :param codimension: The codimension. If codimension = 0, all cofaces + are returned (equivalent of get_star function) + :type codimension: int + :returns: The (simplices of the) cofaces of a simplex + :rtype: list of tuples(simplex, filtration) + """ + cdef vector[int] csimplex = simplex + cdef int num_parameters = self.num_parameters + # for i in simplex: + # csimplex.push_back(i) + cdef vector[simplex_filtration_type] cofaces \ + = self.get_ptr().get_cofaces(csimplex, codimension) + ct = [] + for filtered_simplex in cofaces: + v = [] + for vertex in filtered_simplex.first: + v.append(vertex) + ct.append((v, np.asarray(filtered_simplex.second))) + return ct + + def get_boundaries(self, simplex): + """This function returns a generator with the boundaries of a given N-simplex. + If you do not need the filtration values, the boundary can also be obtained as + :code:`itertools.combinations(simplex,len(simplex)-1)`. + + :param simplex: The N-simplex, represented by a list of vertex. + :type simplex: list of int. + :returns: The (simplices of the) boundary of a simplex + :rtype: generator with tuples(simplex, filtration) + """ + cdef pair[Simplex_tree_multi_boundary_iterator, Simplex_tree_multi_boundary_iterator] it = self.get_ptr().get_boundary_iterators(simplex) + + # while it.first != it.second: + # yield self.get_ptr().get_simplex_and_filtration(dereference(it.first)) + # preincrement(it.first) + cdef int num_parameters = self.get_ptr().get_number_of_parameters() + while it.first != it.second: + # yield self.get_ptr().get_simplex_and_filtration(dereference(it)) + pair = self.get_ptr().get_simplex_and_filtration(dereference(it.first)) + yield (np.asarray(pair.first, dtype=int),np.asarray( pair.second)) + preincrement(it.first) + def remove_maximal_simplex(self, simplex): + """This function removes a given maximal N-simplex from the simplicial + complex. + + :param simplex: The N-simplex, represented by a list of vertex. + :type simplex: list of int + + .. note:: + + The dimension of the simplicial complex may be lower after calling + remove_maximal_simplex than it was before. However, + :func:`upper_bound_dimension` + method will return the old value, which + remains a valid upper bound. If you care, you can call + :func:`dimension` + to recompute the exact dimension. + """ + self.get_ptr().remove_maximal_simplex(simplex) + + # def prune_above_filtration(self, filtration)->bool: + # """Prune above filtration value given as parameter. + + # :param filtration: Maximum threshold value. + # :type filtration: float + # :returns: The filtration modification information. + # :rtype: bool + + + # .. note:: + + # Note that the dimension of the simplicial complex may be lower + # after calling + # :func:`prune_above_filtration` + # than it was before. However, + # :func:`upper_bound_dimension` + # will return the old value, which remains a + # valid upper bound. If you care, you can call + # :func:`dimension` + # method to recompute the exact dimension. + # """ + # return self.get_ptr().prune_above_filtration(filtration) + def prune_above_dimension(self, int dimension): + """Remove all simplices of dimension greater than a given value. + + :param dimension: Maximum dimension value. + :type dimension: int + :returns: The modification information. + :rtype: bool + """ + return self.get_ptr().prune_above_dimension(dimension) + def expansion(self, int max_dim)->SimplexTreeMulti: + """Expands the simplex tree containing only its one skeleton + until dimension max_dim. + + The expanded simplicial complex until dimension :math:`d` + attached to a graph :math:`G` is the maximal simplicial complex of + dimension at most :math:`d` admitting the graph :math:`G` as + :math:`1`-skeleton. + The filtration value assigned to a simplex is the maximal filtration + value of one of its edges. + + The simplex tree must contain no simplex of dimension bigger than + 1 when calling the method. + + :param max_dim: The maximal dimension. + :type max_dim: int + """ + with nogil: + self.get_ptr().expansion(max_dim) + # This is a fix for multipersistence. FIXME expansion in c++ + self.get_ptr().make_filtration_non_decreasing() + return self + + def make_filtration_non_decreasing(self)->bool: + """This function ensures that each simplex has a higher filtration + value than its faces by increasing the filtration values. + + :returns: True if any filtration value was modified, + False if the filtration was already non-decreasing. + :rtype: bool + """ + cdef bool out + with nogil: + out = self.get_ptr().make_filtration_non_decreasing() + return out + + def reset_filtration(self, filtration, min_dim = 0): + """This function resets the filtration value of all the simplices of dimension at least min_dim. Resets all the + simplex tree when `min_dim = 0`. + `reset_filtration` may break the filtration property with `min_dim > 0`, and it is the user's responsibility to + make it a valid filtration (using a large enough `filt_value`, or calling `make_filtration_non_decreasing` + afterwards for instance). + + :param filtration: New threshold value. + :type filtration: float. + :param min_dim: The minimal dimension. Default value is 0. + :type min_dim: int. + """ + self.get_ptr().reset_filtration(Finitely_critical_multi_filtration(filtration), min_dim) + + + + # def extend_filtration(self): + # """ Extend filtration for computing extended persistence. This function only uses the filtration values at the + # 0-dimensional simplices, and computes the extended persistence diagram induced by the lower-star filtration + # computed with these values. + # + # .. note:: + # + # Note that after calling this function, the filtration values are actually modified within the simplex tree. + # The function :func:`extended_persistence` retrieves the original values. + # + # .. note:: + # + # Note that this code creates an extra vertex internally, so you should make sure that the simplex tree does + # not contain a vertex with the largest possible value (i.e., 4294967295). + # + # This `notebook `_ + # explains how to compute an extension of persistence called extended persistence. + # """ + # self.get_ptr().compute_extended_filtration() + + # def extended_persistence(self, homology_coeff_field=11, min_persistence=0): + # """This function retrieves good values for extended persistence, and separate the diagrams into the Ordinary, + # Relative, Extended+ and Extended- subdiagrams. + # + # :param homology_coeff_field: The homology coefficient field. Must be a prime number. Default value is 11. Max is 46337. + # :type homology_coeff_field: int + # :param min_persistence: The minimum persistence value (i.e., the absolute value of the difference between the + # persistence diagram point coordinates) to take into account (strictly greater than min_persistence). + # Default value is 0.0. Sets min_persistence to -1.0 to see all values. + # :type min_persistence: float + # :returns: A list of four persistence diagrams in the format described in :func:`persistence`. The first one is + # Ordinary, the second one is Relative, the third one is Extended+ and the fourth one is Extended-. + # See https://link.springer.com/article/10.1007/s10208-008-9027-z and/or section 2.2 in + # https://link.springer.com/article/10.1007/s10208-017-9370-z for a description of these subtypes. + # + # .. note:: + # + # This function should be called only if :func:`extend_filtration` has been called first! + # + # .. note:: + # + # The coordinates of the persistence diagram points might be a little different than the + # original filtration values due to the internal transformation (scaling to [-2,-1]) that is + # performed on these values during the computation of extended persistence. + # + # This `notebook `_ + # explains how to compute an extension of persistence called extended persistence. + # """ + # cdef vector[pair[int, pair[value_type, value_type]]] persistence_result + # if self.pcohptr != NULL: + # del self.pcohptr + # self.pcohptr = new Simplex_tree_persistence_interface(self.get_ptr(), False) + # self.pcohptr.compute_persistence(homology_coeff_field, -1.) + # return self.pcohptr.compute_extended_persistence_subdiagrams(min_persistence) + + # TODO : cython3 + # def expansion_with_blocker(self, max_dim, blocker_func): + # """Expands the Simplex_tree containing only a graph. Simplices corresponding to cliques in the graph are added + # incrementally, faces before cofaces, unless the simplex has dimension larger than `max_dim` or `blocker_func` + # returns `True` for this simplex. + + # The function identifies a candidate simplex whose faces are all already in the complex, inserts it with a + # filtration value corresponding to the maximum of the filtration values of the faces, then calls `blocker_func` + # with this new simplex (represented as a list of int). If `blocker_func` returns `True`, the simplex is removed, + # otherwise it is kept. The algorithm then proceeds with the next candidate. + + # .. warning:: + # Several candidates of the same dimension may be inserted simultaneously before calling `blocker_func`, so + # if you examine the complex in `blocker_func`, you may hit a few simplices of the same dimension that have + # not been vetted by `blocker_func` yet, or have already been rejected but not yet removed. + + # :param max_dim: Expansion maximal dimension value. + # :type max_dim: int + # :param blocker_func: Blocker oracle. + # :type blocker_func: Callable[[List[int]], bool] + # """ + # self.get_ptr().expansion_with_blockers_callback(max_dim, callback, blocker_func) + + # def persistence(self, homology_coeff_field=11, min_persistence=0, persistence_dim_max = False): + # """This function computes and returns the persistence of the simplicial complex. + # + # :param homology_coeff_field: The homology coefficient field. Must be a + # prime number. Default value is 11. Max is 46337. + # :type homology_coeff_field: int + # :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 + # :param persistence_dim_max: If true, the persistent homology for the + # maximal dimension in the complex is computed. If false, it is + # ignored. Default is false. + # :type persistence_dim_max: bool + # :returns: The persistence of the simplicial complex. + # :rtype: list of pairs(dimension, pair(birth, death)) + # """ + # self.compute_persistence(homology_coeff_field, min_persistence, persistence_dim_max) + # return self.pcohptr.get_persistence() + + +## This function is only meant for the edge collapse interface. + def get_edge_list(self): + cdef edge_list out + with nogil: + out = self.get_ptr().get_edge_list() + return out + + def collapse_edges(self, max_dimension:int=None, num:int=1, progress:bool=False, strong:bool=True, full:bool=False, ignore_warning:bool=False)->SimplexTreeMulti: + """Edge collapse for 1-critical 2-parameter clique complex (see https://arxiv.org/abs/2211.05574). + It uses the code from the github repository https://github.com/aj-alonso/filtration_domination . + + Parameters + ---------- + max_dimension:int + Max simplicial dimension of the complex. Unless specified, keeps the same dimension. + num:int + The number of collapses to do. + strong:bool + Whether to use strong collapses or standard collapses (slower, but may remove more edges) + full:bool + Collapses the maximum number of edges if true, i.e., will do (at most) 100 strong collapses and (at most) 100 non-strong collapses afterward. + progress:bool + If true, shows the progress of the number of collapses. + + WARNING + ------- + - This will destroy all of the k-simplices, with k>=2. Be sure to use this with a clique complex, if you want to preserve the homology >= dimension 1. + - This is for 1 critical simplices, with 2 parameter persistence. + Returns + ------- + self:SimplexTreeMulti + A (smaller) simplex tree that has the same homology over this bifiltration. + + """ + # TODO : find a way to do multiple edge collapses without python conversions. + if num <= 0: + return self + assert self.num_parameters == 2, "Number of parameters has to be 2 to use edge collapses ! This is a limitation of Filtration-domination" + if self.dimension > 1 and not ignore_warning: warn("This method ignores simplices of dimension > 1 !") + + max_dimension = self.dimension if max_dimension is None else max_dimension + + # Retrieves the edge list, and send it to filration_domination + edges = self.get_edge_list() + from gudhi.multiparameter_edge_collapse import _collapse_edge_list + edges = _collapse_edge_list(edges, num=num, full=full, strong=strong, progress=progress) + # Retrieves the collapsed simplicial complex + self._reconstruct_from_edge_list(edges, swap=True, expand_dimension=max_dimension) + return self + def _reconstruct_from_edge_list(self, edges, swap:bool=True, expand_dimension:int=None)->SimplexTreeMulti: + """ + Generates a 1-dimensional copy of self, with the edges given as input. Useful for edge collapses + + Input + ----- + - edges : Iterable[(int,int),(float,float)] ## This is the format of the rust library filtration-domination + - swap : bool + If true, will swap self and the collapsed simplextrees. + - expand_dim : int + expands back the simplextree to this dimension + Ouput + ----- + The reduced SimplexTreeMulti having only these edges. + """ + reduced_tree = SimplexTreeMulti(num_parameters=self.num_parameters) + + ## Adds vertices back, with good filtration + if self.num_vertices > 0: + vertices = np.asarray([splx for splx, f in self.get_skeleton(0)], dtype=int).T + vertices_filtration = np.asarray([f for splx, f in self.get_skeleton(0)], dtype=np.float32) + reduced_tree.insert_batch(vertices, vertices_filtration) + + ## Adds edges again + if self.num_simplices - self.num_vertices > 0: + edges_filtration = np.asarray([f for e,f in edges], dtype=np.float32) + edges = np.asarray([e for e, _ in edges], dtype=int).T + reduced_tree.insert_batch(edges, edges_filtration) + if swap: + # Swaps the simplextrees pointers + self.thisptr, reduced_tree.thisptr = reduced_tree.thisptr, self.thisptr # Swaps self and reduced tree (self is a local variable) + if expand_dimension is not None: + self.expansion(expand_dimension) # Expands back the simplextree to the original dimension. + return self if swap else reduced_tree + + @property + def num_parameters(self)->int: + return self.get_ptr().get_number_of_parameters() + def get_simplices_of_dimension(self, dim:int)->np.ndarray: + return np.asarray(self.get_ptr().get_simplices_of_dimension(dim), dtype=int) + def key(self, simplex:list|np.ndarray): + return self.get_ptr().get_key(simplex) + def set_keys_to_enumerate(self)->None: + self.get_ptr().set_keys_to_enumerate() + return + def set_key(self,simplex:list|np.ndarray, key:int)->None: + self.get_ptr().set_key(simplex, key) + return + + + def to_scc(self, path="scc_dataset.txt", progress:bool=True, overwrite:bool=False, ignore_last_generators:bool=True, strip_comments:bool=False, reverse_block:bool=True, rivet_compatible=False)->None: + """ Create a file with the scc2020 standard, representing the n-filtration of the simplextree. + Link : https://bitbucket.org/mkerber/chain_complex_format/src/master/ + + Parameters + ---------- + path:str + path of the file. + ignore_last_generators:bool = True + If false, will include the filtration values of the last free persistence module. + progress:bool = True + Shows the progress bar. + overwrite:bool = False + If true, will overwrite the previous file if it already exists. + ignore_last_generators:bool=True + If true, does not write the final generators to the file. Rivet ignores them. + reverse_block:bool=True + Some obscure programs reverse the inside-block order. + rivet_compatible:bool=False + Returns a firep (old scc2020) format instead. Only Rivet uses this. + + Returns + ------- + Nothing + """ + ### initialize keys + self.set_keys_to_enumerate() + ### File + from os.path import exists + from os import remove + if exists(path): + if not(overwrite): + raise Exception(f"The file {path} already exists. Use the `overwrite` flag if you want to overwrite.") + remove(path) + file = open(path, "a") + file.write("scc2020\n") if not rivet_compatible else file.write("firep\n") + if not strip_comments and not rivet_compatible: file.write("# Number of parameters\n") + num_parameters = self.get_ptr().get_number_of_parameters() + if rivet_compatible: + assert num_parameters == 2 + file.write("Filtration 1\n") + file.write("Filtration 2\n") + else: + file.write(f"{num_parameters}\n") + if not strip_comments: file.write("# Sizes of generating sets\n") + ## WRITES TSR VARIABLES + tsr:list= [0]*(self.dimension+1) # dimension --- 0 + for splx,f in self.get_simplices(): + dim = len(splx)-1 + tsr[dim] += (int)(len(f) // num_parameters) + tsr.reverse() + file.write(" ".join([str(n) for n in tsr])+"\n") + + ## Adds the boundaries to the dictionnary + tsr + dict_splx_to_firep_number = {} + tsr:list = [[] for _ in range(len(tsr))] # tsr stores simplices vertices, according to dimension, and the dictionnary + for dim in range(self.dimension,-1 , -1): # range(2,-1,-1): + for splx,F in self.get_skeleton(dim): + if len(splx) != dim+1: continue + for b,_ in self.get_boundaries(splx): + if not self.key(b) in dict_splx_to_firep_number: + dict_splx_to_firep_number[self.key(b)] = len(tsr[dim-1]) + tsr[dim-1].append(b) + + ## Adds simplices that are not borders to tsr, i.e., simplices not in the dictionnary + for splx,_ in self.get_simplices(): + if not self.key(splx) in dict_splx_to_firep_number: + tsr[len(splx)-1].append(splx) + ## Writes simplices of tsr to file + dim_range = range(self.dimension,0,-1) if ignore_last_generators else range(self.dimension,-1,-1) + for dim in dim_range: # writes block by block + if not strip_comments: file.write(f"# Block of dimension {dim}\n") + if reverse_block: tsr[dim].reverse() + for splx in tsr[dim]: # for simplices of dimension + F = np.concatenate(self.filtration(splx), axis=0) + nbirth = (int)(len(F)//num_parameters) + for i in range(nbirth): + simplex_filtration = F[i*num_parameters:(i+1)*num_parameters] + file.write(" ".join([str(f) for f in simplex_filtration])) + file.write(" ;") + for b,_ in self.get_boundaries(splx): + file.write(f" {dict_splx_to_firep_number[self.key(b)]}") + file.write("\n") + file.close() + return + + def to_rivet(self, path="rivet_dataset.txt", degree:int|None = None, progress:bool=False, overwrite:bool=False, xbins:int|None=None, ybins:int|None=None)->None: + """ Create a file that can be imported by rivet, representing the filtration of the simplextree. + + Parameters + ---------- + path:str + path of the file. + degree:int + The homological degree to ask rivet to compute. + progress:bool = True + Shows the progress bar. + overwrite:bool = False + If true, will overwrite the previous file if it already exists. + Returns + ------- + Nothing + """ + ... + from os.path import exists + from os import remove + if exists(path): + if not(overwrite): + print(f"The file {path} already exists. Use the `overwrite` flag if you want to overwrite.") + return + remove(path) + file = open(path, "a") + file.write("--datatype bifiltration\n") + file.write(f"--homology {degree}\n") if degree is not None else None + file.write(f"-x {xbins}\n") if xbins is not None else None + file.write(f"-y {ybins}\n") if ybins is not None else None + file.write("--xlabel time of appearance\n") + file.write("--ylabel density\n\n") + from tqdm import tqdm + with tqdm(total=self.num_simplices, position=0, disable = not(progress), desc="Writing simplex to file") as bar: + for dim in range(0,self.dimension+1): # Not sure if dimension sort is necessary for rivet. Check ? + file.write(f"# block of dimension {dim}\n") + for s,F in self.get_skeleton(dim): + if len(s) != dim+1: continue + for i in s: + file.write(str(i) + " ") + file.write("; ") + for f in F: + file.write(str(f) + " ") + file.write("\n") + bar.update(1) + file.close() + return + + + + def _get_filtration_values(self, vector[int] degrees, bool inf_to_nan:bool=False)->Iterable[np.ndarray]: + # cdef vector[int] c_degrees = degrees + cdef intptr_t ptr = self.thisptr + cdef vector[vector[vector[value_type]]] out + with nogil: + out = get_filtration_values(ptr, degrees) + filtrations_values = [np.asarray(filtration) for filtration in out] + # Removes infs + if inf_to_nan: + for i,f in enumerate(filtrations_values): + filtrations_values[i][f == np.inf] = np.nan + filtrations_values[i][f == - np.inf] = np.nan + return filtrations_values + + @staticmethod + def _reduce_grid(filtrations_values,resolutions=None, strategy:str="exact", bool unique=True, some_float _q_factor=1., drop_quantiles=[0,0]): + num_parameters = len(filtrations_values) + if resolutions is None and strategy not in ["exact", "precomputed"]: + raise ValueError("Resolutions must be provided for this strategy.") + elif resolutions is not None: + try: + int(resolutions) + resolutions = [resolutions]*num_parameters + except: + pass + try: + a,b=drop_quantiles + except: + a,b=drop_quantiles,drop_quantiles + + if a != 0 or b != 0: + boxes = np.asarray([np.quantile(filtration, [a, b], axis=1, method='closest_observation') for filtration in filtrations_values]) + min_filtration, max_filtration = np.min(boxes, axis=(0,1)), np.max(boxes, axis=(0,1)) # box, birth/death, filtration + filtrations_values = [ + filtration[(m np.asarray([len(f) for f in F])): + return SimplexTreeMulti._reduce_grid(filtrations_values=filtrations_values, resolutions=resolutions, strategy="quantile",_q_factor=1.5*_q_factor) + elif strategy == "regular": + F = [np.linspace(f.min(),f.max(),num=r) for f,r in zip(filtrations_values, resolutions)] + elif strategy == "regular_closest": + + F = [_todo_regular_closest(f,r, unique) for f,r in zip(filtrations_values, resolutions)] + elif strategy == "partition": + F = [_todo_partition(f,r, unique) for f,r in zip(filtrations_values, resolutions)] + elif strategy == "precomputed": + F=filtrations_values + else: + raise Exception("Invalid strategy. Pick either regular, regular_closest, partition, quantile, precomputed or exact.") + + return F + + def get_filtration_grid(self, resolution:Iterable[int]|None=None, degrees:Iterable[int]|None=None, drop_quantiles:float|tuple=0, grid_strategy:str="exact")->Iterable[np.ndarray]: + """ + Returns a grid over the n-filtration, from the simplextree. Usefull for grid_squeeze. TODO : multicritical + + Parameters + ---------- + resolution: list[int] + resolution of the grid, for each parameter + box=None : pair[list[float]] + Grid bounds. format : [low bound, high bound] + If None is given, will use the filtration bounds of the simplextree. + grid_strategy="regular" : string + Either "regular", "quantile", or "exact". + Returns + ------- + List of filtration values, for each parameter, defining the grid. + """ + if degrees is None: + degrees = range(self.dimension+1) + + + ## preprocesses the filtration values: + filtrations_values = np.concatenate(self._get_filtration_values(degrees, inf_to_nan=True), axis=1) + # removes duplicate + sort (nan at the end) + filtrations_values = [np.unique(filtration) for filtration in filtrations_values] + # removes nan + filtrations_values = [filtration[:-1] if np.isnan(filtration[-1]) else filtration for filtration in filtrations_values] + + return self._reduce_grid(filtrations_values=filtrations_values, resolutions=resolution,strategy=grid_strategy,drop_quantiles=drop_quantiles) + + + + def grid_squeeze(self, filtration_grid:np.ndarray|list|None=None, bool coordinate_values=True, force=False, **filtration_grid_kwargs)->SimplexTreeMulti: + """ + Fit the filtration of the simplextree to a grid. + + :param filtration_grid: The grid on which to squeeze. An example of grid can be given by the `get_filtration_grid` method. + :type filtration_grid: list[list[float]] + :param coordinate_values: If true, the filtrations values of the simplices will be set to the coordinate of the filtration grid. + :type coordinate_values: bool + """ + if not force and self._is_squeezed: + raise Exception("SimplexTree already squeezed, use `force=True` if that's really what you want to do.") + #TODO : multi-critical + if filtration_grid is None: filtration_grid = self.get_filtration_grid(**filtration_grid_kwargs) + cdef vector[vector[value_type]] c_filtration_grid = filtration_grid + cdef intptr_t ptr = self.thisptr + if coordinate_values: + self.filtration_grid = c_filtration_grid + with nogil: + squeeze_filtration(ptr, c_filtration_grid, coordinate_values) + return self + + @property + def _is_squeezed(self)->bool: + return self.num_vertices > 0 and self.filtration_grid[0].size() > 0 + + def filtration_bounds(self, degrees:Iterable[int]|None=None, q:float|tuple=0, split_dimension:bool=False)->np.ndarray: + """ + Returns the filtrations bounds of the finite filtration values. + """ + try: + a,b =q + except: + a,b,=q,q + degrees = range(self.dimension+1) if degrees is None else degrees + filtrations_values = self._get_filtration_values(degrees, inf_to_nan=True) ## degree, parameter, pt + boxes = np.array([np.nanquantile(filtration, [a, 1-b], axis=1) for filtration in filtrations_values],dtype=float) + if split_dimension: return boxes + return np.asarray([np.nanmin(boxes, axis=(0,1)), np.nanmax(boxes, axis=(0,1))]) # box, birth/death, filtration + + + + + def fill_lowerstar(self, F, parameter:int)->SimplexTreeMulti: + """ Fills the `dimension`th filtration by the lower-star filtration defined by F. + + Parameters + ---------- + F:1d array + The density over the vertices, that induces a lowerstar filtration. + parameter:int + Which filtration parameter to fill. /!\ python starts at 0. + + Returns + ------- + self:SimplexTreeMulti + """ + # for s, sf in self.get_simplices(): + # self.assign_filtration(s, [f if i != dimension else np.max(np.array(F)[s]) for i,f in enumerate(sf)]) + cdef int c_parameter = parameter + cdef vector[value_type] c_F = F + with nogil: + self.get_ptr().fill_lowerstar(c_F, c_parameter) + return self + + + def project_on_line(self, parameter:int=0, basepoint:None|list|np.ndarray= None)->SimplexTree: + """Converts an multi simplextree to a gudhi simplextree. + Parameters + ---------- + parameter:int = 0 + The parameter to keep. WARNING will crash if the multi simplextree is not well filled. + basepoint:None + Instead of keeping a single parameter, will consider the filtration defined by the diagonal line crossing the basepoint. + WARNING + ------- + There are no safeguard yet, it WILL crash if asking for a parameter that is not filled. + Returns + ------- + A SimplexTree with chosen 1D filtration. + """ + # FIXME : deal with multicritical filtrations + import gudhi as gd + new_simplextree = gd.SimplexTree() + assert parameter < self.get_ptr().get_number_of_parameters() + cdef int c_parameter = parameter + cdef intptr_t old_ptr = self.thisptr + cdef intptr_t new_ptr = new_simplextree.thisptr + cdef vector[value_type] c_basepoint = [] if basepoint is None else basepoint + if basepoint is None: + with nogil: + flatten_from_ptr(old_ptr, new_ptr, c_parameter) + else: + with nogil: + flatten_diag_from_ptr(old_ptr, new_ptr, c_basepoint, c_parameter) + return new_simplextree + + def linear_projections(self, linear_forms:np.ndarray)->Iterable[SimplexTree]: + """ + Compute the 1-parameter projections, w.r.t. given the linear forms, of this simplextree. + + Input + ----- + - Array of shape (num_linear_forms, num_parameters) + + Output + ------ + - List of projected (gudhi) simplextrees. + """ + cdef Py_ssize_t num_projections = linear_forms.shape[0] + cdef Py_ssize_t num_parameters = linear_forms.shape[1] + if num_projections == 0: return [] + cdef vector[vector[value_type]] c_linear_forms = linear_forms + assert num_parameters==self.num_parameters, f"The linear forms has to have the same number of parameter as the simplextree ({self.num_parameters})." + + # Gudhi copies are faster than inserting simplices one by one + import gudhi as gd + flattened_simplextree = gd.SimplexTree() + cdef intptr_t multi_prt = self.thisptr + cdef intptr_t flattened_ptr = flattened_simplextree.thisptr + with nogil: + flatten_from_ptr(multi_prt, flattened_ptr, num_parameters) + out = [flattened_simplextree] + [gd.SimplexTree(flattened_simplextree) for _ in range(num_projections-1)] + + # Fills the 1-parameter simplextrees. + cdef vector[intptr_t] out_ptrs = [st.thisptr for st in out] + with nogil: + for i in range(num_projections): + linear_projection_from_ptr(out_ptrs[i], multi_prt, c_linear_forms[i]) + return out + + + def set_num_parameter(self, num:int): + """ + Sets the numbers of parameters. + WARNING : it will resize all the filtrations to this size. + """ + self.get_ptr().resize_all_filtrations(num) + self.get_ptr().set_number_of_parameters(num) + return + + def __eq__(self, other:SimplexTreeMulti): + """Test for structural equality + :returns: True if the 2 simplex trees are equal, False otherwise. + :rtype: bool + """ + return dereference(self.get_ptr()) == dereference(other.get_ptr()) + +cdef intptr_t _get_copy_intptr(SimplexTreeMulti stree) nogil: + return (new Simplex_tree_multi_interface(dereference(stree.get_ptr()))) + + +def _todo_regular_closest(cnp.ndarray[some_float,ndim=1] f, int r, bool unique): + f_regular = np.linspace(np.min(f),np.max(f),num=r) + f_regular_closest = np.asarray([f[np.argmin(np.abs(f-x))] for x in f_regular]) + if unique: f_regular_closest = np.unique(f_regular_closest) + return f_regular_closest + +def _todo_partition(cnp.ndarray[some_float,ndim=1] data,int resolution, bool unique): + if data.shape[0] < resolution: resolution=data.shape[0] + k = data.shape[0] // resolution + partitions = np.partition(data, k) + f = partitions[[i*k for i in range(resolution)]] + if unique: f= np.unique(f) + return f + + + +def _simplextree_multify(simplextree:SimplexTree, num_parameters:int=2, default_values=[])->SimplexTreeMulti: + """Converts a gudhi simplextree to a multi simplextree. + Parameters + ---------- + parameters:int = 2 + The number of filtrations + Returns + ------- + A multi simplextree, with first filtration value being the one from the original simplextree. + """ + if isinstance(simplextree, SimplexTreeMulti): + return simplextree + st = SimplexTreeMulti(num_parameters=num_parameters) + cdef int c_num_parameters = num_parameters + cdef intptr_t old_ptr = simplextree.thisptr + cdef intptr_t new_ptr = st.thisptr + cdef vector[value_type] c_default_values=default_values + with nogil: + multify_from_ptr(old_ptr, new_ptr, c_num_parameters, c_default_values) + return st + +def _safe_simplextree_multify(simplextree:SimplexTree,int num_parameters=2)->SimplexTreeMulti: + if isinstance(simplextree, SimplexTreeMulti): + return simplextree + simplices = [[] for _ in range(simplextree.dimension()+1)] + filtration_values = [[] for _ in range(simplextree.dimension()+1)] + for s,f in simplextree.get_simplices(): + filtration_values[len(s)-1].append(f) + simplices[len(s)-1].append(s) + st_multi = SimplexTreeMulti(num_parameters=1) + for batch_simplices, batch_filtrations in zip(simplices,filtration_values): + st_multi.insert_batch(np.asarray(batch_simplices).T, np.asarray(batch_filtrations)[:,None]) + if num_parameters > 1: + st_multi.set_num_parameter(num_parameters) + return st_multi \ No newline at end of file diff --git a/src/python/include/Simplex_tree_interface_multi.h b/src/python/include/Simplex_tree_interface_multi.h new file mode 100644 index 0000000000..9c0e14bcf0 --- /dev/null +++ b/src/python/include/Simplex_tree_interface_multi.h @@ -0,0 +1,254 @@ +/* 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): Vincent Rouvreau + * + * Copyright (C) 2016 Inria + * + * Modification(s): + * - 2022/11 David Loiseaux, Hannah Schreiber : adapt for multipersistence. + * - YYYY/MM Author: Description of the modification + */ + +#pragma once + + +#include "Simplex_tree_interface.h" +#include "Simplex_tree_multi.h" +#include "multi_filtrations/finitely_critical_filtrations.h" + +#include +#include +#include // std::pair +#include +#include // for std::distance +#include + +namespace Gudhi::multiparameter { + +template +class Simplex_tree_interface_multi : public Simplex_tree_interface { + public: + using Python_filtration_type = std::vector; // TODO : std::conditional + using Base = Simplex_tree; + using Filtration_value = typename Base::Filtration_value; + using Vertex_handle = typename Base::Vertex_handle; + using Simplex_handle = typename Base::Simplex_handle; + using Insertion_result = typename std::pair; + using Simplex = std::vector; + using Simplex_and_filtration = std::pair; + using Filtered_simplices = std::vector; + 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; + typedef bool (*blocker_func_t)(Simplex simplex, void *user_data); + using euler_chars_type = std::vector; + + public: + + Extended_filtration_data efd; + + bool find_simplex(const Simplex& simplex) { + return (Base::find(simplex) != Base::null_simplex()); + } + + void assign_simplex_filtration(const Simplex& simplex, const Filtration_value& filtration) { + Base::assign_filtration(Base::find(simplex), filtration); + Base::clear_filtration(); + } + + + bool insert(const Simplex& simplex, const Filtration_value& filtration ) { + Insertion_result result = Base::insert_simplex_and_subfaces(simplex, filtration); + if (result.first != Base::null_simplex()) + Base::clear_filtration(); + return (result.second); + } + + // Do not interface this function, only used in alpha complex interface for complex creation + bool insert_simplex(const Simplex& simplex, const Filtration_value& filtration ) { + Insertion_result result = Base::insert_simplex(simplex, filtration); + return (result.second); + } +// bool insert_simplex(const Simplex& simplex, const Python_filtration_type& filtration ) { +// Filtration_value& filtration_ = *(Filtration_value*)(&filtration); // Jardinage for no copy. +// Insertion_result result = Base::insert_simplex(simplex, filtration); +// return (result.second); +// } + + // Do not interface this function, only used in interface for complex creation + bool insert_simplex_and_subfaces(const Simplex& simplex, const Filtration_value& filtration ) { + Insertion_result result = Base::insert_simplex_and_subfaces(simplex, filtration); + return (result.second); + } +// bool insert_simplex_and_subfaces(const Simplex& simplex, const Python_filtration_type& filtration ) { +// Filtration_value& filtration_ = *(Filtration_value*)(&filtration); // Jardinage for no copy. +// Insertion_result result = Base::insert_simplex_and_subfaces(simplex, filtration); +// return (result.second); +// } + + // Do not interface this function, only used in strong witness interface for complex creation + bool insert_simplex(const std::vector& simplex, const Filtration_value& filtration ) { + Insertion_result result = Base::insert_simplex(simplex, filtration); + return (result.second); + } + + // Do not interface this function, only used in strong witness interface for complex creation + bool insert_simplex_and_subfaces(const std::vector& simplex, const Filtration_value& filtration) { + Insertion_result result = Base::insert_simplex_and_subfaces(simplex, filtration); + return (result.second); + } + typename SimplexTreeOptions::value_type* simplex_filtration(const Simplex& simplex) { + auto& filtration = Base::filtration_mutable(Base::find(simplex)); + return &filtration[0]; // We return the pointer to get a numpy view afterward + } + + + 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); // why not push back ? + // } + auto it = Base::simplex_vertex_range(f_simplex); + Simplex simplex(it.begin(), it.end()); + std::reverse(simplex.begin(), simplex.end()); + return std::make_pair(std::move(simplex), &Base::filtration_mutable(f_simplex)[0]); + } + + Filtered_simplices get_star(const Simplex& simplex) { + Filtered_simplices star; + for (auto f_simplex : Base::star_simplex_range(Base::find(simplex))) { + Simplex simplex_star; + for (auto vertex : Base::simplex_vertex_range(f_simplex)) { + simplex_star.insert(simplex_star.begin(), vertex); + } + star.push_back(std::make_pair(simplex_star, &Base::filtration_mutable(f_simplex)[0])); + } + return star; + } + + Filtered_simplices get_cofaces(const Simplex& simplex, int dimension) { + Filtered_simplices cofaces; + for (auto f_simplex : Base::cofaces_simplex_range(Base::find(simplex), dimension)) { + Simplex simplex_coface; + for (auto vertex : Base::simplex_vertex_range(f_simplex)) { + simplex_coface.insert(simplex_coface.begin(), vertex); + } + cofaces.push_back(std::make_pair(simplex_coface, &Base::filtration_mutable(f_simplex)[0])); + } + return cofaces; + } + + void compute_extended_filtration() { + throw std::logic_error("Incompatible with multipers"); + } + + Simplex_tree_interface_multi* collapse_edges(int nb_collapse_iteration) { + throw std::logic_error("Incompatible with multipers"); + } + + + +// ######################## MULTIPERS STUFF + void set_keys_to_enumerate(){ + int count = 0; + for (auto sh : Base::filtration_simplex_range()) + Base::assign_key(sh, count++); + } + + int get_key(const Simplex& simplex){ + return Base::key(Base::find(simplex)); + } + + void set_key(const Simplex& simplex, int key){ + Base::assign_key(Base::find(simplex), key); + return; + } + + // Fills a parameter with a lower-star filtration + void fill_lowerstar(const std::vector& filtration, int axis){ + using value_type=options_multi::value_type; + for (auto &SimplexHandle : Base::complex_simplex_range()){ + std::vector& current_birth = Base::filtration_mutable(SimplexHandle); + value_type to_assign = -1*std::numeric_limits::infinity(); + for (auto vertex : Base::simplex_vertex_range(SimplexHandle)){ + to_assign = std::max(filtration[vertex], to_assign); + } + current_birth[axis] = to_assign; + // Base::assign_filtration(SimplexHandle, current_birth); + } + } + + + using simplices_list = std::vector>; + simplices_list get_simplices_of_dimension(int dimension){ + simplices_list simplex_list; + simplex_list.reserve(Base::num_simplices()); + for (auto simplexhandle : Base::skeleton_simplex_range(dimension)){ + if (Base::dimension(simplexhandle) == dimension){ + std::vector simplex; + simplex.reserve(dimension+1); + for (int vertex : Base::simplex_vertex_range(simplexhandle)) + simplex.push_back(vertex); + simplex_list.push_back(simplex); + } + } +/* simplex_list.shrink_to_fit();*/ + return simplex_list; + } + using edge_list = std::vector, std::pair>>; + edge_list get_edge_list(){ + edge_list simplex_list; + simplex_list.reserve(Base::num_simplices()); + for (auto &simplexHandle : Base::skeleton_simplex_range(1)){ + if (Base::dimension(simplexHandle) == 1){ + std::pair simplex; + auto it = Base::simplex_vertex_range(simplexHandle).begin(); + simplex = {*it, *(++it)}; + const auto& f = Base::filtration(simplexHandle); + simplex_list.push_back({simplex, {f[0], f[1]}}); + } + } +/* simplex_list.shrink_to_fit();*/ + return simplex_list; + } + void resize_all_filtrations(int num){ //TODO : that is for 1 critical filtrations + if (num < 0) return; + for(const auto &SimplexHandle : Base::complex_simplex_range()){ + Base::filtration_mutable(SimplexHandle).resize(num);; + } + } + + +}; + + +using interface_std = Simplex_tree; // Interface not necessary (smaller so should do less segfaults) +using interface_multi = Simplex_tree_interface_multi; + +// Wrappers of the functions in Simplex_tree_multi.h, to deal with the "pointer only" python interface +void flatten_diag_from_ptr(const uintptr_t splxptr, const uintptr_t newsplxptr, const std::vector basepoint, int dimension){ // for python + auto &st = get_simplextree_from_pointer(newsplxptr); + auto &st_multi = get_simplextree_from_pointer(splxptr); + flatten_diag(st,st_multi,basepoint, dimension); +} +void multify_from_ptr(uintptr_t splxptr,uintptr_t newsplxptr, const int dimension, const multi_filtration_type& default_values){ //for python + auto &st = get_simplextree_from_pointer(splxptr); + auto &st_multi = get_simplextree_from_pointer(newsplxptr); + multify(st, st_multi, dimension, default_values); +} +void flatten_from_ptr(uintptr_t splxptr, uintptr_t newsplxptr, const int dimension = 0){ // for python + auto &st = get_simplextree_from_pointer(newsplxptr); + auto &st_multi = get_simplextree_from_pointer(splxptr); + flatten(st, st_multi, dimension); +} +template +void linear_projection_from_ptr(const uintptr_t ptr, const uintptr_t ptr_multi, Args...args){ + auto &st = get_simplextree_from_pointer(ptr); + auto &st_multi = get_simplextree_from_pointer(ptr_multi); + linear_projection(st, st_multi, args...); +} + + +} // namespace Gudhi + diff --git a/src/python/include/Simplex_tree_multi.h b/src/python/include/Simplex_tree_multi.h new file mode 100644 index 0000000000..2725346e5e --- /dev/null +++ b/src/python/include/Simplex_tree_multi.h @@ -0,0 +1,238 @@ +/* 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): Hannah Schreiber + * + * Copyright (C) 2014 Inria """ + * + * + * Modification(s): + * - 2022/11 David Loiseaux / Hannah Schreiber : added multify / flatten to interface standard simplextree. + * - YYYY/MM Author: Description of the modification + */ +#ifndef SIMPLEX_TREE_MULTI_H_ +#define SIMPLEX_TREE_MULTI_H_ + +#include +#include +#include "multi_filtrations/finitely_critical_filtrations.h" +#include "multi_filtrations/line.h" + + + + + +namespace Gudhi::multiparameter { +/** Model of SimplexTreeOptions. + * + * Maximum number of simplices to compute persistence is std::numeric_limits::max() + * (about 4 billions of simplices). */ + + +struct Simplex_tree_options_multidimensional_filtration { +public: + typedef linear_indexing_tag Indexing_tag; + typedef int Vertex_handle; + typedef float value_type; + using Filtration_value = multi_filtrations::Finitely_critical_multi_filtration; + typedef std::uint32_t Simplex_key; + static const bool store_key = true; + static const bool store_filtration = true; + static const bool contiguous_vertices = false; + static const bool link_nodes_by_label = true; + static const bool stable_simplex_handles = false; + static const bool is_multi_parameter = true; +}; + + + +using options_multi = Simplex_tree_options_multidimensional_filtration; +using options_std = Simplex_tree_options_full_featured; +using simplextree_std = Simplex_tree; +using simplextree_multi = Simplex_tree; + +using multi_filtration_type = std::vector; +using multi_filtration_grid = std::vector; + +// Retrieves a simplextree from a pointer. As the python simplex_tree and simplex_tree_multi don't know each other we have to exchange them using pointers. +template +simplextreeinterface& get_simplextree_from_pointer(const uintptr_t splxptr){ //DANGER + simplextreeinterface &st = *(simplextreeinterface*)(splxptr); + return st; +} + +// Turns a 1-parameter simplextree into a multiparameter simplextree, and keeps the 1-filtration in the 1st axis +template +void multify(simplextree_std &st, simplextree_multi &st_multi, const int num_parameters, const typename simplextree_multi::Options::Filtration_value& default_values={}){ + typename simplextree_multi::Options::Filtration_value f(num_parameters); + for (auto i = 0u; i(default_values.size()), static_cast(num_parameters-1));i++) + f[i+1] = default_values[i]; + std::vector simplex; + simplex.reserve(st.dimension()+1); + for (auto &simplex_handle : st.complex_simplex_range()){ + simplex.clear(); + for (auto vertex : st.simplex_vertex_range(simplex_handle)) + simplex.push_back(vertex); + + if (num_parameters > 0) + f[0] = st.filtration(simplex_handle); + auto filtration_copy = f; + st_multi.insert_simplex(simplex,filtration_copy); + + } +} + + + +// Turns a multi-parameter simplextree into a 1-parameter simplextree +template +void flatten(simplextree_std &st, simplextree_multi &st_multi, const int dimension = 0){ + for (const auto &simplex_handle : st_multi.complex_simplex_range()){ + std::vector simplex; + for (auto vertex : st_multi.simplex_vertex_range(simplex_handle)) + simplex.push_back(vertex); + typename simplextree_multi::Options::value_type f = dimension >= 0 ? st_multi.filtration(simplex_handle)[dimension] : 0; + st.insert_simplex(simplex,f); + } +} + +// Applies a linear form (i.e. scalar product with Riesz rpz) to the filtration to flatten a simplextree multi +template +void linear_projection(simplextree_std &st, simplextree_multi &st_multi, const std::vector& linear_form){ + auto sh = st.complex_simplex_range().begin(); + auto sh_multi = st_multi.complex_simplex_range().begin(); + auto end = st.complex_simplex_range().end(); + typename simplextree_multi::Options::Filtration_value multi_filtration; + for (; sh != end; ++sh, ++sh_multi){ + multi_filtration = st_multi.filtration(*sh_multi); + auto projected_filtration = multi_filtration.linear_projection(linear_form); + st.assign_filtration(*sh, projected_filtration); + } +} + +template +void flatten_diag(simplextree_std &st, simplextree_multi &st_multi, const std::vector basepoint, int dimension){ + multi_filtrations::Line l(basepoint); + for (const auto &simplex_handle : st_multi.complex_simplex_range()){ + std::vector simplex; + for (auto vertex : st_multi.simplex_vertex_range(simplex_handle)) + simplex.push_back(vertex); + + std::vector f = st_multi.filtration(simplex_handle); + if (dimension <0) dimension = 0; + typename simplextree_multi::Options::value_type new_filtration = l.push_forward(f)[dimension]; + st.insert_simplex(simplex,new_filtration); + } +} + + + + + +/// @brief turns filtration value x into coordinates in the grid +/// @tparam out_type +/// @param x +/// @param grid +/// @return +template +inline void find_coordinates(vector_like& x, const multi_filtration_grid &grid){ + // TODO: optimize with, e.g., dichotomy + + for (auto parameter = 0u; parameter < grid.size(); parameter++){ + const auto& filtration = grid[parameter]; // assumes its sorted + const auto to_project = x[parameter]; + if constexpr (std::numeric_limits::has_infinity) + if (to_project == std::numeric_limits::infinity()){ + x[parameter] = std::numeric_limits::infinity(); + continue; + } + if (to_project >= filtration.back()){ + x[parameter] = filtration.size()-1; + continue; + } // deals with infinite value at the end of the grid + + unsigned int i = 0; + while (to_project > filtration[i] && i &st_multi = *(Gudhi::Simplex_tree*)(splxptr); + auto num_parameters = static_cast(st_multi.get_number_of_parameters()); + if (grid.size() != num_parameters){ + std::cerr << "Bad grid !" << std::endl; + throw; + } + for (const auto &simplex_handle : st_multi.complex_simplex_range()){ + auto& simplex_filtration = st_multi.filtration_mutable(simplex_handle); + find_coordinates(simplex_filtration, grid); // turns the simplexfiltration into coords in the grid + if (!coordinate_values){ + for (auto parameter = 0u; parameter < num_parameters; parameter++) + simplex_filtration[parameter] = grid[parameter][simplex_filtration[parameter]]; + } + } + return; +} + +// retrieves the filtration values of a simplextree. Useful to generate a grid. +std::vector get_filtration_values(const uintptr_t splxptr, const std::vector °rees){ + Simplex_tree &st_multi = *(Gudhi::Simplex_tree*)(splxptr); + int num_parameters = st_multi.get_number_of_parameters(); + std::vector out(degrees.size(), multi_filtration_grid(num_parameters)); + std::vector degree_index(degrees.size()); + int count = 0; + for (auto degree : degrees){ + degree_index[degree] = count; count++; + out[degree_index[degree]].reserve(st_multi.num_simplices()); + } + + for (const auto &simplex_handle : st_multi.complex_simplex_range()){ + const auto filtration = st_multi.filtration(simplex_handle); + const auto degree = st_multi.dimension(simplex_handle); + if (std::find(degrees.begin(), degrees.end(), degree) == degrees.end()) continue; + for (int parameter=0; parameter < num_parameters; parameter++){ + out[degree_index[degree]][parameter].push_back(filtration[parameter]); + } + } + return out; + +} + +} // namespace Gudhi + + + +namespace std { + +template<> +class numeric_limits +{ +public: + static Gudhi::multiparameter::multi_filtration_type infinity() throw(){ + return Gudhi::multiparameter::multi_filtration_type(1, std::numeric_limits::infinity()); + }; + + + static Gudhi::multiparameter::multi_filtration_type quiet_NaN() throw(){ + return Gudhi::multiparameter::multi_filtration_type(1, numeric_limits::quiet_NaN()); + }; + +}; + +} + + +#endif // SIMPLEX_TREE_MULTI_H_ diff --git a/src/python/include/multi_filtrations/box.h b/src/python/include/multi_filtrations/box.h new file mode 100644 index 0000000000..a4c28fdf1f --- /dev/null +++ b/src/python/include/multi_filtrations/box.h @@ -0,0 +1,164 @@ +/* This file is part of the MMA Library - https://gitlab.inria.fr/dloiseau/multipers - which is released under MIT. + * See file LICENSE for full license details. + * Author(s): David Loiseaux + * + * Copyright (C) 2021 Inria + * + * Modification(s): + * + */ +/** + * @file box.h + * @author David Loiseaux + * @brief BOX. + */ + +#ifndef BOX_H_INCLUDED +#define BOX_H_INCLUDED + +#include +#include +#include +#include +#include + +#include "finitely_critical_filtrations.h" + + + + +/** + * @brief Simple box in $\mathbb R^n$ . + */ + +namespace Gudhi::multiparameter::multi_filtrations{ + +template +class Box +{ + + using point_type = Finitely_critical_multi_filtration; +public: + Box(); + Box(T a, T b, T c, T d) : bottomCorner_({a,b}), upperCorner_({c,d}) {}; + Box(const point_type& bottomCorner, const point_type& upperCorner); + Box(const std::pair& box); + + void inflate(T delta); + const point_type& get_bottom_corner() const; + const point_type& get_upper_corner() const; + point_type& get_bottom_corner(); + point_type& get_upper_corner(); + bool contains(const point_type& point) const; + void infer_from_filters(const std::vector &Filters_list); + bool is_trivial() const ; + std::pair get_pair() const{ + return {bottomCorner_,upperCorner_}; + } + std::pair get_pair(){ + return {bottomCorner_,upperCorner_}; + } + +private: + point_type bottomCorner_; + point_type upperCorner_; +}; + +template +inline Box::Box() +{} + +template +inline Box::Box(const point_type &bottomCorner, const point_type &upperCorner) + : bottomCorner_(bottomCorner), + upperCorner_(upperCorner) +{ + assert(bottomCorner.size() == upperCorner.size() + && bottomCorner <= upperCorner + && "This box is trivial !"); +} + +template +inline Box::Box(const std::pair &box) + : bottomCorner_(box.first), + upperCorner_(box.second) +{} + + +template +inline void Box::inflate(T delta) +{ + bottomCorner_ -= delta; + upperCorner_ += delta; +} + +template +inline void Box::infer_from_filters(const std::vector &Filters_list){ + int dimension = Filters_list.size(); + int nsplx = Filters_list[0].size(); + std::vector lower(dimension); + std::vector upper(dimension); + for (int i = 0; i < dimension; i++){ + T min = Filters_list[i][0]; + T max = Filters_list[i][0]; + for (int j=1; j +inline bool Box::is_trivial() const { + return bottomCorner_.empty() || upperCorner_.empty() || bottomCorner_.size() != upperCorner_.size(); +} + +template +inline const typename Box::point_type &Box::get_bottom_corner() const +{ + return bottomCorner_; +} + +template +inline const typename Box::point_type &Box::get_upper_corner() const +{ + return upperCorner_; +} + +template +inline typename Box::point_type &Box::get_bottom_corner() +{ + return bottomCorner_; +} + +template +inline typename Box::point_type &Box::get_upper_corner() +{ + return upperCorner_; +} + +template +inline bool Box::contains(const point_type &point) const +{ + if (point.size() != bottomCorner_.size()) return false; + + return bottomCorner_ <= point && point <= upperCorner_; +} + +template +std::ostream& operator<<(std::ostream& os, const Box& box) +{ + os << "Box -- Bottom corner : "; + os << box.get_bottom_corner(); + os << ", Top corner : "; + os << box.get_upper_corner(); + return os; +} + +} // namespace Gudhi + + +#endif // BOX_H_INCLUDED diff --git a/src/python/include/multi_filtrations/finitely_critical_filtrations.h b/src/python/include/multi_filtrations/finitely_critical_filtrations.h new file mode 100644 index 0000000000..a75c45c641 --- /dev/null +++ b/src/python/include/multi_filtrations/finitely_critical_filtrations.h @@ -0,0 +1,150 @@ +#ifndef FINITELY_CRITICAL_FILTRATIONS_H_ +#define FINITELY_CRITICAL_FILTRATIONS_H_ + +#include +#include +#include + +namespace Gudhi::multiparameter::multi_filtrations{ + +template +class Finitely_critical_multi_filtration : public std::vector { + // Class to prevent doing illegal stuff with the standard library, e.g., compare two vectors +public: + Finitely_critical_multi_filtration() : std::vector() {}; + Finitely_critical_multi_filtration(int n) : std::vector(n, -std::numeric_limits::infinity()) {}; // minus infinity by default + Finitely_critical_multi_filtration(int n, T value) : std::vector(n,value) {}; + Finitely_critical_multi_filtration(std::initializer_list init) : std::vector(init) {}; + Finitely_critical_multi_filtration(const std::vector& v) : std::vector(v) {}; + Finitely_critical_multi_filtration(typename std::vector::iterator it_begin,typename std::vector::iterator it_end) : std::vector(it_begin, it_end) {}; + Finitely_critical_multi_filtration(typename std::vector::const_iterator it_begin,typename std::vector::const_iterator it_end) : std::vector(it_begin, it_end) {}; + + + operator std::vector&() const { + return *this; + } + std::vector get_vector() const{ + return static_cast>(*this); + } + + + + friend bool operator<(const Finitely_critical_multi_filtration& a, const Finitely_critical_multi_filtration& b) + { + bool isSame = true; + int n = std::min(a.size(), b.size()); + for (int i = 0; i < n; ++i){ + if (a[i] > b[i]) return false; + if (isSame && a[i] != b[i]) isSame = false; + } + if (isSame) return false; + return true; + } + friend bool operator<=(const Finitely_critical_multi_filtration& a, const Finitely_critical_multi_filtration& b) + { + int n = std::min(a.size(), b.size()); + for (int i = 0; i < n; ++i){ + if (a[i] > b[i]) return false; + } + return true; + } + + + + //GREATER THAN OPERATORS + friend bool operator>(const Finitely_critical_multi_filtration& a, const Finitely_critical_multi_filtration& b) + { + return b=(const Finitely_critical_multi_filtration& a, const Finitely_critical_multi_filtration& b) + { + return b<=a; + } + + Finitely_critical_multi_filtration& operator=(const Finitely_critical_multi_filtration& a){ + std::vector::operator=(a); + return *this; + } + + std::vector& _convert_back(){ + return *this; + } + + + + + friend Finitely_critical_multi_filtration& operator-=(Finitely_critical_multi_filtration &result, const Finitely_critical_multi_filtration &to_substract){ + std::transform(result.begin(), result.end(), to_substract.begin(),result.begin(), std::minus()); + return result; + } + friend Finitely_critical_multi_filtration& operator+=(Finitely_critical_multi_filtration &result, const Finitely_critical_multi_filtration &to_add){ + std::transform(result.begin(), result.end(), to_add.begin(),result.begin(), std::plus()); + return result; + } + + friend Finitely_critical_multi_filtration& operator-=(Finitely_critical_multi_filtration &result, const T &to_substract){ + // std::transform(result.begin(), result.end(), to_substract.begin(),result.begin(), std::minus()); + for (auto & truc : result){ + truc -= to_substract; + } + return result; + } + friend Finitely_critical_multi_filtration& operator+=(Finitely_critical_multi_filtration &result, const T &to_add){ + for (auto & truc : result){ + truc += to_add; + } + return result; + } + + static std::vector> to_python(const std::vector>& to_convert){ + return std::vector>(to_convert.begin(), to_convert.end()); + } + + + static std::vector> from_python(const std::vector>& to_convert){ + return std::vector>(to_convert.begin(), to_convert.end());; + } + void push_to(const Finitely_critical_multi_filtration& x){ + if (this->size() != x.size()){ + std::cerr << "Does only work with 1-critical filtrations ! Sizes " << this->size() << " and " << x.size() << "are different !" << std::endl; + throw std::logic_error("Bad sizes"); + } + for (unsigned int i = 0; i < x.size(); i++) + this->at(i) = this->at(i) > x[i] ? this->at(i) : x[i]; + } + // Warning, this function assumes that the comparisons checks have already been made ! + void insert_new(Finitely_critical_multi_filtration to_concatenate){ + this->insert( + this->end(), std::move_iterator(to_concatenate.begin()), std::move_iterator(to_concatenate.end()) + ); + } + + + // scalar product of a filtration value with x. + T linear_projection(const std::vector& x){ + T projection=0; + unsigned int size = std::min(x.size(), this->size()); + for (auto i =0u; iat(i); + return projection; + } + + // easy debug + friend std::ostream& operator<<(std::ostream& stream, const Finitely_critical_multi_filtration& truc){ + stream << "["; + for(unsigned int i = 0; i < truc.size()-1; i++){ + stream << truc[i] << ", "; + } + if(!truc.empty()) stream << truc.back(); + stream << "]"; + return stream; + } + + + + +}; + + +} // namespace Gudhi +#endif // FINITELY_CRITICAL_FILTRATIONS_H_ diff --git a/src/python/include/multi_filtrations/line.h b/src/python/include/multi_filtrations/line.h new file mode 100644 index 0000000000..0fa82cf715 --- /dev/null +++ b/src/python/include/multi_filtrations/line.h @@ -0,0 +1,94 @@ +/* This file is part of the MMA Library - https://gitlab.inria.fr/dloiseau/multipers - which is released under MIT. + * See file LICENSE for full license details. + * Author(s): David Loiseaux + * + * Copyright (C) 2022 Inria + * + * Modification(s): + */ +/** + * @file line_filtration_translation.h + * @author David Loiseaux + * @brief + */ + + +#ifndef LINE_FILTRATION_TRANSLATION_H_INCLUDED +#define LINE_FILTRATION_TRANSLATION_H_INCLUDED + +#include "box.h" +#include "finitely_critical_filtrations.h" + +namespace Gudhi::multiparameter::multi_filtrations{ + + + template + class Line + { + + public: + using point_type = Finitely_critical_multi_filtration; + Line(); + Line(point_type x); + Line(point_type x, point_type v); + point_type push_forward(point_type x) const; + point_type push_back(point_type x) const; + int get_dim() const; + std::pair get_bounds(const Box &box) const; + + + private: + point_type basepoint_; // any point on the line + point_type direction_; // direction of the line + + }; + template + Line::Line(){} + + template + Line::Line(point_type x){ + this->basepoint_.swap(x); + // this->direction_ = {}; // diagonal line + } + template + Line::Line(point_type x, point_type v){ + this->basepoint_.swap(x); + this->direction_.swap(v); + } + template + typename Line::point_type Line::push_forward(point_type x) const { //TODO remove copy + x -= basepoint_; + T t = - std::numeric_limits::infinity();; + for (std::size_t i = 0; idirection_.size() > i ? direction_[i] : 1; + t = std::max(t, x[i]/dir); + } + point_type out(basepoint_.size()); + for (unsigned int i = 0; i < out.size(); i++) + out[i] = basepoint_[i] + t * (this->direction_.size() > i ? direction_[i] : 1) ; + return out; + } + template + typename Line::point_type Line::push_back(point_type x) const{ + x -= basepoint_; + T t = std::numeric_limits::infinity(); + for (unsigned int i = 0; idirection_.size() > i ? direction_[i] : 1; + t = std::min(t, x[i]/dir); + } + point_type out(basepoint_.size()); + for (unsigned int i = 0; i < out.size(); i++) + out[i] = basepoint_[i] + t * (this->direction_.size() > i ? direction_[i] : 1) ; + return out; + } + template + int Line::get_dim() const{ + return basepoint_.size(); + } + template + std::pair::point_type, typename Line::point_type> Line::get_bounds(const Box &box) const{ + return {this->push_forward(box.get_bottom_corner()), this->push_back(box.get_upper_corner())}; + } +} + +#endif // LINE_FILTRATION_TRANSLATION_H_INCLUDED