diff options
Diffstat (limited to 'src/Simplex_tree/include/gudhi/Simplex_tree.h')
| -rw-r--r-- | src/Simplex_tree/include/gudhi/Simplex_tree.h | 431 |
1 files changed, 354 insertions, 77 deletions
diff --git a/src/Simplex_tree/include/gudhi/Simplex_tree.h b/src/Simplex_tree/include/gudhi/Simplex_tree.h index fafdb01c..4177a0b8 100644 --- a/src/Simplex_tree/include/gudhi/Simplex_tree.h +++ b/src/Simplex_tree/include/gudhi/Simplex_tree.h @@ -24,6 +24,9 @@ #include <boost/iterator/transform_iterator.hpp> #include <boost/graph/adjacency_list.hpp> #include <boost/range/adaptor/reversed.hpp> +#include <boost/range/adaptor/transformed.hpp> +#include <boost/range/size.hpp> +#include <boost/container/static_vector.hpp> #ifdef GUDHI_USE_TBB #include <tbb/parallel_sort.h> @@ -41,6 +44,24 @@ namespace Gudhi { +/** \addtogroup simplex_tree + * @{ + */ + +/** + * \class Extended_simplex_type Simplex_tree.h gudhi/Simplex_tree.h + * \brief Extended simplex type data structure for representing the type of simplices in an extended filtration. + * + * \details The extended simplex type can be either UP (which means + * that the simplex was present originally, and is thus part of the ascending extended filtration), DOWN (which means + * that the simplex is the cone of an original simplex, and is thus part of the descending extended filtration) or + * EXTRA (which means the simplex is the cone point). + * + * Details may be found in \cite Cohen-Steiner2009 and section 2.2 in \cite Carriere16. + * + */ +enum class Extended_simplex_type {UP, DOWN, EXTRA}; + struct Simplex_tree_options_full_featured; /** @@ -82,10 +103,11 @@ class Simplex_tree { // Simplex_key next to each other). typedef typename boost::container::flat_map<Vertex_handle, Node> Dictionary; - /* \brief Set of nodes sharing a same parent in the simplex tree. */ - /* \brief Set of nodes sharing a same parent in the simplex tree. */ + /** \brief Set of nodes sharing a same parent in the simplex tree. */ typedef Simplex_tree_siblings<Simplex_tree, Dictionary> Siblings; + + struct Key_simplex_base_real { Key_simplex_base_real() : key_(-1) {} void assign_key(Simplex_key k) { key_ = k; } @@ -99,6 +121,12 @@ class Simplex_tree { void assign_key(Simplex_key); Simplex_key key() const; }; + struct Extended_filtration_data { + Filtration_value minval; + Filtration_value maxval; + Extended_filtration_data(){} + Extended_filtration_data(Filtration_value vmin, Filtration_value vmax): minval(vmin), maxval(vmax) {} + }; typedef typename std::conditional<Options::store_key, Key_simplex_base_real, Key_simplex_base_dummy>::type Key_simplex_base; @@ -119,7 +147,10 @@ class Simplex_tree { public: /** \brief Handle type to a simplex contained in the simplicial complex represented - * by the simplex tree. */ + * by the simplex tree. + * + * They are essentially pointers into internal vectors, and any insertion or removal + * of a simplex may invalidate any other Simplex_handle in the complex. */ typedef typename Dictionary::iterator Simplex_handle; private: @@ -161,6 +192,12 @@ class Simplex_tree { typedef Simplex_tree_boundary_simplex_iterator<Simplex_tree> Boundary_simplex_iterator; /** \brief Range over the simplices of the boundary of a simplex. */ typedef boost::iterator_range<Boundary_simplex_iterator> Boundary_simplex_range; + /** \brief Iterator over the simplices of the boundary of a simplex and their opposite vertices. + * + * 'value_type' is std::pair<Simplex_handle, Vertex_handle>. */ + typedef Simplex_tree_boundary_opposite_vertex_simplex_iterator<Simplex_tree> Boundary_opposite_vertex_simplex_iterator; + /** \brief Range over the simplices of the boundary of a simplex and their opposite vertices. */ + typedef boost::iterator_range<Boundary_opposite_vertex_simplex_iterator> Boundary_opposite_vertex_simplex_range; /** \brief Iterator over the simplices of the simplicial complex. * * 'value_type' is Simplex_handle. */ @@ -232,11 +269,9 @@ class Simplex_tree { * * The filtration must be valid. If the filtration has not been initialized yet, the * method initializes it (i.e. order the simplices). If the complex has changed since the last time the filtration - * was initialized, please call `initialize_filtration()` to recompute it. */ + * was initialized, please call `clear_filtration()` or `initialize_filtration()` to recompute it. */ Filtration_simplex_range const& filtration_simplex_range(Indexing_tag = Indexing_tag()) { - if (filtration_vect_.empty()) { - initialize_filtration(); - } + maybe_initialize_filtration(); return filtration_vect_; } @@ -246,8 +281,8 @@ class Simplex_tree { * which is consequenlty * equal to \f$(-1)^{\text{dim} \sigma}\f$ the canonical orientation on the simplex. */ - Simplex_vertex_range simplex_vertex_range(Simplex_handle sh) { - assert(sh != null_simplex()); // Empty simplex + Simplex_vertex_range simplex_vertex_range(Simplex_handle sh) const { + GUDHI_CHECK(sh != null_simplex(), "empty simplex"); return Simplex_vertex_range(Simplex_vertex_iterator(this, sh), Simplex_vertex_iterator(this)); } @@ -272,6 +307,23 @@ class Simplex_tree { Boundary_simplex_iterator(this)); } + /** \brief Given a simplex, returns a range over the simplices of its boundary and their opposite vertices. + * + * The boundary of a simplex is the set of codimension \f$1\f$ subsimplices of the simplex. + * If the simplex is \f$[v_0, \cdots ,v_d]\f$, with canonical orientation induced by \f$ v_0 < \cdots < v_d \f$, the + * iterator enumerates the simplices of the boundary in the order: + * \f$[v_0,\cdots,\widehat{v_i},\cdots,v_d]\f$ for \f$i\f$ from \f$d\f$ to \f$0\f$, where \f$\widehat{v_i}\f$ means + * that the vertex \f$v_i\f$, known as the opposite vertex, is omitted from boundary, but returned as the second + * element of a pair. + * + * @param[in] sh Simplex for which the boundary is computed. + */ + template<class SimplexHandle> + Boundary_opposite_vertex_simplex_range boundary_opposite_vertex_simplex_range(SimplexHandle sh) { + return Boundary_opposite_vertex_simplex_range(Boundary_opposite_vertex_simplex_iterator(this, sh), + Boundary_opposite_vertex_simplex_iterator(this)); + } + /** @} */ // end range and iterator methods /** \name Constructor/Destructor * @{ */ @@ -286,7 +338,7 @@ class Simplex_tree { /** \brief User-defined copy constructor reproduces the whole tree structure. */ Simplex_tree(const Simplex_tree& complex_source) { #ifdef DEBUG_TRACES - std::cout << "Simplex_tree copy constructor" << std::endl; + std::clog << "Simplex_tree copy constructor" << std::endl; #endif // DEBUG_TRACES copy_from(complex_source); } @@ -296,7 +348,7 @@ class Simplex_tree { */ Simplex_tree(Simplex_tree && complex_source) { #ifdef DEBUG_TRACES - std::cout << "Simplex_tree move constructor" << std::endl; + std::clog << "Simplex_tree move constructor" << std::endl; #endif // DEBUG_TRACES move_from(complex_source); @@ -313,7 +365,7 @@ class Simplex_tree { /** \brief User-defined copy assignment reproduces the whole tree structure. */ Simplex_tree& operator= (const Simplex_tree& complex_source) { #ifdef DEBUG_TRACES - std::cout << "Simplex_tree copy assignment" << std::endl; + std::clog << "Simplex_tree copy assignment" << std::endl; #endif // DEBUG_TRACES // Self-assignment detection if (&complex_source != this) { @@ -330,7 +382,7 @@ class Simplex_tree { */ Simplex_tree& operator=(Simplex_tree&& complex_source) { #ifdef DEBUG_TRACES - std::cout << "Simplex_tree move assignment" << std::endl; + std::clog << "Simplex_tree move assignment" << std::endl; #endif // DEBUG_TRACES // Self-assignment detection if (&complex_source != this) { @@ -450,10 +502,19 @@ class Simplex_tree { return true; } + /** \brief Returns the filtration value of a simplex. + * + * Same as `filtration()`, but does not handle `null_simplex()`. + */ + static Filtration_value filtration_(Simplex_handle sh) { + GUDHI_CHECK (sh != null_simplex(), "null simplex"); + return sh->second.filtration(); + } + public: /** \brief Returns the key associated to a simplex. * - * The filtration must be initialized. + * If no key has been assigned, returns `null_key()`. * \pre SimplexTreeOptions::store_key */ static Simplex_key key(Simplex_handle sh) { @@ -463,7 +524,6 @@ class Simplex_tree { /** \brief Returns the simplex that has index idx in the filtration. * * The filtration must be initialized. - * \pre SimplexTreeOptions::store_key */ Simplex_handle simplex(Simplex_key idx) const { return filtration_vect_[idx]; @@ -499,8 +559,7 @@ class Simplex_tree { return Dictionary_it(nullptr); } - /** \brief Returns a key different for all keys associated to the - * simplices of the simplicial complex. */ + /** \brief Returns a fixed number not in the interval [0, `num_simplices()`). */ static Simplex_key null_key() { return -1; } @@ -645,10 +704,10 @@ class Simplex_tree { return true; } - private: - /** \brief Inserts a simplex represented by a vector of vertex. - * @param[in] simplex vector of Vertex_handles, representing the vertices of the new simplex. The vector must be - * sorted by increasing vertex handle order. + protected: + /** \brief Inserts a simplex represented by a range of vertex. + * @param[in] simplex range of Vertex_handles, representing the vertices of the new simplex. The range must be + * sorted by increasing vertex handle order, and not empty. * @param[in] filtration the filtration value assigned to the new simplex. * @return If the new simplex is inserted successfully (i.e. it was not in the * simplicial complex yet) the bool is set to true and the Simplex_handle is the handle assigned @@ -660,12 +719,13 @@ class Simplex_tree { * null_simplex. * */ - std::pair<Simplex_handle, bool> insert_vertex_vector(const std::vector<Vertex_handle>& simplex, + template <class RandomVertexHandleRange = std::initializer_list<Vertex_handle>> + std::pair<Simplex_handle, bool> insert_simplex_raw(const RandomVertexHandleRange& simplex, Filtration_value filtration) { Siblings * curr_sib = &root_; std::pair<Simplex_handle, bool> res_insert; auto vi = simplex.begin(); - for (; vi != simplex.end() - 1; ++vi) { + for (; vi != std::prev(simplex.end()); ++vi) { GUDHI_CHECK(*vi != null_vertex(), "cannot use the dummy null_vertex() as a real vertex"); res_insert = curr_sib->members_.emplace(*vi, Node(curr_sib, filtration)); if (!(has_children(res_insert.first))) { @@ -686,9 +746,10 @@ class Simplex_tree { return std::pair<Simplex_handle, bool>(null_simplex(), false); } // otherwise the insertion has succeeded - size is a size_type - if (static_cast<int>(simplex.size()) - 1 > dimension_) { + int dim = static_cast<int>(boost::size(simplex)) - 1; + if (dim > dimension_) { // Update dimension if needed - dimension_ = static_cast<int>(simplex.size()) - 1; + dimension_ = dim; } return res_insert; } @@ -729,7 +790,7 @@ class Simplex_tree { // Copy before sorting std::vector<Vertex_handle> copy(first, last); std::sort(std::begin(copy), std::end(copy)); - return insert_vertex_vector(copy, filtration); + return insert_simplex_raw(copy, filtration); } /** \brief Insert a N-simplex and all his subfaces, from a N-simplex represented by a range of @@ -755,12 +816,7 @@ class Simplex_tree { if (first == last) return { null_simplex(), true }; // FIXME: false would make more sense to me. - // Copy before sorting - // Thread local is not available on XCode version < V.8 - It will slow down computation -#ifdef GUDHI_CAN_USE_CXX11_THREAD_LOCAL - thread_local -#endif // GUDHI_CAN_USE_CXX11_THREAD_LOCAL - std::vector<Vertex_handle> copy; + thread_local std::vector<Vertex_handle> copy; copy.clear(); copy.insert(copy.end(), first, last); std::sort(copy.begin(), copy.end()); @@ -827,7 +883,7 @@ class Simplex_tree { /** Returns the Siblings containing a simplex.*/ template<class SimplexHandle> - Siblings* self_siblings(SimplexHandle sh) { + static Siblings* self_siblings(SimplexHandle sh) { if (sh->second.children()->parent() == sh->first) return sh->second.children()->oncles(); else @@ -850,15 +906,13 @@ class Simplex_tree { } public: - /** \brief Initializes the filtrations, i.e. sort the - * simplices according to their order in the filtration and initializes all Simplex_keys. + /** \brief Initializes the filtration cache, i.e. sorts the + * simplices according to their order in the filtration. * - * After calling this method, filtration_simplex_range() becomes valid, and each simplex is - * assigned a Simplex_key corresponding to its order in the filtration (from 0 to m-1 for a - * simplicial complex with m simplices). + * It always recomputes the cache, even if one already exists. * - * Will be automatically called when calling filtration_simplex_range() - * if the filtration has never been initialized yet. */ + * Any insertion, deletion or change of filtration value invalidates this cache, + * which can be cleared with clear_filtration(). */ void initialize_filtration() { filtration_vect_.clear(); filtration_vect_.reserve(num_simplices()); @@ -880,6 +934,21 @@ class Simplex_tree { std::stable_sort(filtration_vect_.begin(), filtration_vect_.end(), is_before_in_filtration(this)); #endif } + /** \brief Initializes the filtration cache if it isn't initialized yet. + * + * Automatically called by filtration_simplex_range(). */ + void maybe_initialize_filtration() { + if (filtration_vect_.empty()) { + initialize_filtration(); + } + } + /** \brief Clears the filtration cache produced by initialize_filtration(). + * + * Useful when initialize_filtration() has already been called and we perform an operation + * (say an insertion) that invalidates the cache. */ + void clear_filtration() { + filtration_vect_.clear(); + } private: /** Recursive search of cofaces @@ -903,7 +972,7 @@ class Simplex_tree { // If we reached the end of the vertices, and the simplex has more vertices than the given simplex // => we found a coface - // Add a coface if we wan't the star or if the number of vertices of the current simplex matches with nbVertices + // Add a coface if we want the star or if the number of vertices of the current simplex matches with nbVertices bool addCoface = (star || curr_nbVertices == nbVertices); if (addCoface) cofaces.push_back(simplex); @@ -1021,8 +1090,8 @@ class Simplex_tree { * * Inserts all vertices and edges given by a OneSkeletonGraph. * OneSkeletonGraph must be a model of - * <a href="http://www.boost.org/doc/libs/1_65_1/libs/graph/doc/EdgeListGraph.html">boost::EdgeListGraph</a> - * and <a href="http://www.boost.org/doc/libs/1_65_1/libs/graph/doc/PropertyGraph.html">boost::PropertyGraph</a>. + * <a href="https://www.boost.org/doc/libs/release/libs/graph/doc/VertexAndEdgeListGraph.html">boost::VertexAndEdgeListGraph</a> + * and <a href="https://www.boost.org/doc/libs/release/libs/graph/doc/PropertyGraph.html">boost::PropertyGraph</a>. * * The vertex filtration value is accessible through the property tag * vertex_filtration_t. @@ -1042,7 +1111,10 @@ class Simplex_tree { // the simplex tree must be empty assert(num_simplices() == 0); - if (boost::num_vertices(skel_graph) == 0) { + // is there a better way to let the compiler know that we don't mean Simplex_tree::num_vertices? + using boost::num_vertices; + + if (num_vertices(skel_graph) == 0) { return; } if (num_edges(skel_graph) == 0) { @@ -1051,25 +1123,21 @@ class Simplex_tree { dimension_ = 1; } - root_.members_.reserve(boost::num_vertices(skel_graph)); + root_.members_.reserve(num_vertices(skel_graph)); // probably useless in most cases + auto verts = vertices(skel_graph) | boost::adaptors::transformed([&](auto v){ + return Dit_value_t(v, Node(&root_, get(vertex_filtration_t(), skel_graph, v))); }); + root_.members_.insert(boost::begin(verts), boost::end(verts)); + // This automatically sorts the vertices, the graph concept doesn't guarantee the order in which we iterate. - typename boost::graph_traits<OneSkeletonGraph>::vertex_iterator v_it, - v_it_end; - for (std::tie(v_it, v_it_end) = boost::vertices(skel_graph); v_it != v_it_end; - ++v_it) { - root_.members_.emplace_hint( - root_.members_.end(), *v_it, - Node(&root_, boost::get(vertex_filtration_t(), skel_graph, *v_it))); - } std::pair<typename boost::graph_traits<OneSkeletonGraph>::edge_iterator, - typename boost::graph_traits<OneSkeletonGraph>::edge_iterator> boost_edges = boost::edges(skel_graph); + typename boost::graph_traits<OneSkeletonGraph>::edge_iterator> boost_edges = edges(skel_graph); // boost_edges.first is the equivalent to boost_edges.begin() // boost_edges.second is the equivalent to boost_edges.end() for (; boost_edges.first != boost_edges.second; boost_edges.first++) { auto edge = *(boost_edges.first); auto u = source(edge, skel_graph); auto v = target(edge, skel_graph); - if (u == v) throw "Self-loops are not simplicial"; + if (u == v) throw std::invalid_argument("Self-loops are not simplicial"); // We cannot skip edges with the wrong orientation and expect them to // come a second time with the right orientation, that does not always // happen in practice. emplace() should be a NOP when an element with the @@ -1084,10 +1152,25 @@ class Simplex_tree { } sh->second.children()->members().emplace(v, - Node(sh->second.children(), boost::get(edge_filtration_t(), skel_graph, edge))); + Node(sh->second.children(), get(edge_filtration_t(), skel_graph, edge))); } } + /** \brief Inserts several vertices. + * @param[in] vertices A range of Vertex_handle + * @param[in] filt filtration value of the new vertices (the same for all) + * + * This may be faster than inserting the vertices one by one, especially in a random order. + * 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 <class VertexRange> + void insert_batch_vertices(VertexRange const& vertices, Filtration_value filt = 0) { + 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)); + if (dimension_ < 0 && !root_.members_.empty()) dimension_ = 0; + } + /** \brief Expands the Simplex_tree containing only its one skeleton * until dimension max_dim. * @@ -1101,6 +1184,7 @@ class Simplex_tree { * 1 when calling the method. */ void expansion(int max_dim) { if (max_dim <= 1) return; + clear_filtration(); // Drop the cache. dimension_ = max_dim; for (Dictionary_it root_it = root_.members_.begin(); root_it != root_.members_.end(); ++root_it) { @@ -1123,10 +1207,7 @@ class Simplex_tree { Dictionary_it next = siblings->members().begin(); ++next; -#ifdef GUDHI_CAN_USE_CXX11_THREAD_LOCAL - thread_local -#endif // GUDHI_CAN_USE_CXX11_THREAD_LOCAL - std::vector<std::pair<Vertex_handle, Node> > inter; + thread_local std::vector<std::pair<Vertex_handle, Node> > inter; for (Dictionary_it s_h = siblings->members().begin(); s_h != siblings->members().end(); ++s_h, ++next) { Simplex_handle root_sh = find_vertex(s_h->first); @@ -1243,6 +1324,7 @@ class Simplex_tree { Siblings * new_sib = new Siblings(siblings, // oncles simplex->first, // parent boost::adaptors::reverse(intersection)); // boost::container::ordered_unique_range_t + simplex->second.assign_children(new_sib); std::vector<Vertex_handle> blocked_new_sib_vertex_list; // As all intersections are inserted, we can call the blocker function on all new_sib members for (auto new_sib_member = new_sib->members().begin(); @@ -1265,7 +1347,6 @@ class Simplex_tree { new_sib->members().erase(blocked_new_sib_member); } // ensure recursive call - simplex->second.assign_children(new_sib); siblings_expansion_with_blockers(new_sib, max_dim, k - 1, block_simplex); } } else { @@ -1275,7 +1356,7 @@ class Simplex_tree { } } - /* \private Returns the Simplex_handle composed of the vertex list (from the Simplex_handle), plus the given + /** \private Returns the Simplex_handle composed of the vertex list (from the Simplex_handle), plus the given * Vertex_handle if the Vertex_handle is found in the Simplex_handle children list. * Returns null_simplex() if it does not exist */ @@ -1314,9 +1395,8 @@ class Simplex_tree { /** \brief This function ensures that each simplex has a higher filtration value than its faces by increasing the * filtration values. * @return True if any filtration value was modified, false if the filtration was already non-decreasing. - * \post Some simplex tree functions require the filtration to be valid. `make_filtration_non_decreasing()` - * function is not launching `initialize_filtration()` but returns the filtration modification information. If the - * complex has changed , please call `initialize_filtration()` to recompute it. + * + * If a simplex has a `NaN` filtration value, it is considered lower than any other defined filtration value. */ bool make_filtration_non_decreasing() { bool modified = false; @@ -1326,6 +1406,8 @@ class Simplex_tree { modified |= rec_make_filtration_non_decreasing(simplex.second.children()); } } + if(modified) + clear_filtration(); // Drop the cache. return modified; } @@ -1347,7 +1429,9 @@ class Simplex_tree { }); Filtration_value max_filt_border_value = filtration(*max_border); - if (simplex.second.filtration() < max_filt_border_value) { + // Replacing if(f<max) with 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 (!(simplex.second.filtration() >= max_filt_border_value)) { // Store the filtration modification information modified = true; simplex.second.assign_filtration(max_filt_border_value); @@ -1363,22 +1447,22 @@ class Simplex_tree { public: /** \brief Prune above filtration value given as parameter. * @param[in] filtration Maximum threshold value. - * @return The filtration modification information. - * \post Some simplex tree functions require the filtration to be valid. `prune_above_filtration()` - * function is not launching `initialize_filtration()` but returns the filtration modification information. If the - * complex has changed , please call `initialize_filtration()` to recompute it. + * @return True if any simplex was removed, false if all simplices already had a value below the threshold. * \post Note that the dimension of the simplicial complex may be lower after calling `prune_above_filtration()` * 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) { - return rec_prune_above_filtration(root(), filtration); + bool modified = rec_prune_above_filtration(root(), filtration); + if(modified) + clear_filtration(); // Drop the cache. + return modified; } private: bool rec_prune_above_filtration(Siblings* sib, Filtration_value filt) { auto&& list = sib->members(); - auto last = std::remove_if(list.begin(), list.end(), [=](Dit_value_t& simplex) { + auto last = std::remove_if(list.begin(), list.end(), [this,filt](Dit_value_t& simplex) { if (simplex.second.filtration() <= filt) return false; if (has_children(&simplex)) rec_delete(simplex.second.children()); // dimension may need to be lowered @@ -1418,14 +1502,14 @@ class Simplex_tree { for (Simplex_handle sh : complex_simplex_range()) { #ifdef DEBUG_TRACES for (auto vertex : simplex_vertex_range(sh)) { - std::cout << " " << vertex; + std::clog << " " << vertex; } - std::cout << std::endl; + std::clog << std::endl; #endif // DEBUG_TRACES int sh_dimension = dimension(sh); if (sh_dimension >= dimension_) - // Stop browsing as soon as the dimension is reached, no need to go furter + // Stop browsing as soon as the dimension is reached, no need to go further return false; new_dimension = (std::max)(new_dimension, sh_dimension); } @@ -1439,7 +1523,6 @@ class Simplex_tree { * @param[in] sh Simplex handle on the maximal simplex to remove. * \pre Please check the simplex has no coface before removing it. * \exception std::invalid_argument In debug mode, if sh has children. - * \post Be aware that removing is shifting data in a flat_map (initialize_filtration to be done). * \post Note that the dimension of the simplicial complex may be lower after calling `remove_maximal_simplex()` * 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. @@ -1465,6 +1548,200 @@ class Simplex_tree { } } + /** \brief Retrieve the original filtration value for a given simplex in the Simplex_tree. Since the + * computation of extended persistence requires modifying the filtration values, this function can be used + * to recover the original values. Moreover, computing extended persistence requires adding new simplices + * in the Simplex_tree. Hence, this function also outputs the type of each simplex. It can be either UP (which means + * that the simplex was present originally, and is thus part of the ascending extended filtration), DOWN (which means + * that the simplex is the cone of an original simplex, and is thus part of the descending extended filtration) or + * EXTRA (which means the simplex is the cone point). See the definition of Extended_simplex_type. Note that if the simplex type is DOWN, the original filtration value + * is set to be the original filtration value of the corresponding (not coned) original simplex. + * \pre This function should be called only if `extend_filtration()` has been called first! + * \post The output filtration value is supposed to be the same, but might be a little different, than the + * original filtration value, due to the internal transformation (scaling to [-2,-1]) that is + * performed on the original filtration values during the computation of extended persistence. + * @param[in] f Filtration value of the simplex in the extended (i.e., modified) filtration. + * @param[in] efd Structure containing the minimum and maximum values of the original filtration. This the output of `extend_filtration()`. + * @return A pair containing the original filtration value of the simplex as well as the simplex type. + */ + std::pair<Filtration_value, Extended_simplex_type> decode_extended_filtration(Filtration_value f, const Extended_filtration_data& efd){ + std::pair<Filtration_value, Extended_simplex_type> p; + Filtration_value minval = efd.minval; + Filtration_value maxval = efd.maxval; + if (f >= -2 && f <= -1){ + p.first = minval + (maxval-minval)*(f + 2); p.second = Extended_simplex_type::UP; + } + else if (f >= 1 && f <= 2){ + p.first = minval - (maxval-minval)*(f - 2); p.second = Extended_simplex_type::DOWN; + } + else{ + p.first = std::numeric_limits<Filtration_value>::quiet_NaN(); p.second = Extended_simplex_type::EXTRA; + } + return p; + }; + + /** \brief 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. + * \post Note that after calling this function, the filtration + * values are actually modified. The function `decode_extended_filtration()` + * retrieves the original values and outputs the extended simplex type. + * \pre 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 Vertex_handle. + * @return A data structure containing the maximum and minimum values of the original filtration. + * It is meant to be provided as input to `decode_extended_filtration()` in order to retrieve + * the original filtration values for each simplex. + */ + Extended_filtration_data extend_filtration() { + clear_filtration(); // Drop the cache. + + // Compute maximum and minimum of filtration values + Vertex_handle maxvert = std::numeric_limits<Vertex_handle>::min(); + Filtration_value minval = std::numeric_limits<Filtration_value>::infinity(); + Filtration_value maxval = -std::numeric_limits<Filtration_value>::infinity(); + for (auto sh = root_.members().begin(); sh != root_.members().end(); ++sh){ + Filtration_value f = this->filtration(sh); + minval = std::min(minval, f); + maxval = std::max(maxval, f); + maxvert = std::max(sh->first, maxvert); + } + + GUDHI_CHECK(maxvert < std::numeric_limits<Vertex_handle>::max(), std::invalid_argument("Simplex_tree contains a vertex with the largest Vertex_handle")); + maxvert += 1; + + Simplex_tree st_copy = *this; + + // Add point for coning the simplicial complex + this->insert_simplex_raw({maxvert}, -3); + + // For each simplex + std::vector<Vertex_handle> vr; + for (auto sh_copy : st_copy.complex_simplex_range()){ + + // Locate simplex + vr.clear(); + for (auto vh : st_copy.simplex_vertex_range(sh_copy)){ + vr.push_back(vh); + } + auto sh = this->find(vr); + + // Create cone on simplex + vr.push_back(maxvert); + if (this->dimension(sh) == 0){ + Filtration_value v = this->filtration(sh); + Filtration_value scaled_v = (v-minval)/(maxval-minval); + // Assign ascending value between -2 and -1 to vertex + this->assign_filtration(sh, -2 + scaled_v); + // Assign descending value between 1 and 2 to cone on vertex + this->insert_simplex(vr, 2 - scaled_v); + } + else{ + // Assign value -3 to simplex and cone on simplex + this->assign_filtration(sh, -3); + this->insert_simplex(vr, -3); + } + } + + // Automatically assign good values for simplices + this->make_filtration_non_decreasing(); + + // Return the filtration data + Extended_filtration_data efd(minval, maxval); + return efd; + } + + /** \brief Returns a vertex of `sh` that has the same filtration value as `sh` if it exists, and `null_vertex()` otherwise. + * + * For a lower-star filtration built with `make_filtration_non_decreasing()`, this is a way to invert the process and find out which vertex had its filtration value propagated to `sh`. + * If several vertices have the same filtration value, the one it returns is arbitrary. */ + Vertex_handle vertex_with_same_filtration(Simplex_handle sh) { + auto filt = filtration_(sh); + for(auto v : simplex_vertex_range(sh)) + if(filtration_(find_vertex(v)) == filt) + return v; + return null_vertex(); + } + + /** \brief Returns an edge of `sh` that has the same filtration value as `sh` if it exists, and `null_simplex()` otherwise. + * + * For a flag-complex built with `expansion()`, this is a way to invert the process and find out which edge had its filtration value propagated to `sh`. + * If several edges have the same filtration value, the one it returns is arbitrary. + * + * \pre `sh` must have dimension at least 1. */ + Simplex_handle edge_with_same_filtration(Simplex_handle sh) { + // See issue #251 for potential speed improvements. + auto&& vertices = simplex_vertex_range(sh); // vertices in decreasing order + auto end = std::end(vertices); + auto vi = std::begin(vertices); + GUDHI_CHECK(vi != end, "empty simplex"); + auto v0 = *vi; + ++vi; + GUDHI_CHECK(vi != end, "simplex of dimension 0"); + if(std::next(vi) == end) return sh; // shortcut for dimension 1 + boost::container::static_vector<Vertex_handle, 40> suffix; + suffix.push_back(v0); + auto filt = filtration_(sh); + do + { + Vertex_handle v = *vi; + auto&& children1 = find_vertex(v)->second.children()->members_; + for(auto w : suffix){ + // Can we take advantage of the fact that suffix is ordered? + Simplex_handle s = children1.find(w); + if(filtration_(s) == filt) + return s; + } + suffix.push_back(v); + } + while(++vi != end); + return null_simplex(); + } + + /** \brief Returns a minimal face of `sh` that has the same filtration value as `sh`. + * + * For a filtration built with `make_filtration_non_decreasing()`, this is a way to invert the process and find out which simplex had its filtration value propagated to `sh`. + * If several minimal (for inclusion) simplices have the same filtration value, the one it returns is arbitrary, and it is not guaranteed to be the one with smallest dimension. */ + Simplex_handle minimal_simplex_with_same_filtration(Simplex_handle sh) { + auto filt = filtration_(sh); + // Naive implementation, it can be sped up. + for(auto b : boundary_simplex_range(sh)) + if(filtration_(b) == filt) + return minimal_simplex_with_same_filtration(b); + return sh; // None of its faces has the same filtration. + } + + public: + /** \brief 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[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) { + rec_reset_filtration(&root_, filt_value, min_dim); + clear_filtration(); // Drop the cache. + } + + private: + /** \brief Recursively resets filtration value when minimal depth <= 0. + * @param[in] sib Siblings to be parsed. + * @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) { + for (auto sh = sib->members().begin(); sh != sib->members().end(); ++sh) { + if (min_depth <= 0) { + sh->second.assign_filtration(filt_value); + } + if (has_children(sh)) { + rec_reset_filtration(sh->second.children(), filt_value, min_depth - 1); + } + } + } + private: Vertex_handle null_vertex_; /** \brief Total number of simplices in the complex, without the empty simplex.*/ @@ -1542,7 +1819,7 @@ struct Simplex_tree_options_fast_persistence { static const bool contiguous_vertices = true; }; -/** @} */ // end defgroup simplex_tree +/** @}*/ // end addtogroup simplex_tree } // namespace Gudhi |