diff options
Diffstat (limited to 'src/Simplex_tree/include/gudhi/Simplex_tree.h')
| -rw-r--r-- | src/Simplex_tree/include/gudhi/Simplex_tree.h | 810 |
1 files changed, 810 insertions, 0 deletions
diff --git a/src/Simplex_tree/include/gudhi/Simplex_tree.h b/src/Simplex_tree/include/gudhi/Simplex_tree.h new file mode 100644 index 00000000..5caa5e25 --- /dev/null +++ b/src/Simplex_tree/include/gudhi/Simplex_tree.h @@ -0,0 +1,810 @@ + /* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * Author(s): Clément Maria + * + * Copyright (C) 2014 INRIA Sophia Antipolis-Méditerranée (France) + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + */ + +#ifndef GUDHI_SIMPLEX_TREE_H +#define GUDHI_SIMPLEX_TREE_H + +#include <algorithm> +#include <boost/container/flat_map.hpp> +#include <boost/iterator/transform_iterator.hpp> +#include <boost/graph/adjacency_list.hpp> +#include "gudhi/Simplex_tree/Simplex_tree_node_explicit_storage.h" +#include "gudhi/Simplex_tree/Simplex_tree_siblings.h" +#include "gudhi/Simplex_tree/Simplex_tree_iterators.h" +#include "gudhi/Simplex_tree/indexing_tag.h" + +namespace Gudhi{ + +/** \defgroup simplex_tree Filtered Complexes Package + * + * A simplicial complex \f$\mathbf{K}\f$ + * on a set of vertices \f$V = \{1, \cdots ,|V|\}\f$ is a collection of simplices + * \f$\{\sigma\}\f$, + * \f$\sigma \subseteq V\f$ such that \f$\tau \subseteq \sigma \in \mathbf{K} \rightarrow \tau \in + * \mathbf{K}\f$. The + * dimension \f$n=|\sigma|-1\f$ of \f$\sigma\f$ is its number of elements minus \f$1\f$. + * + * A filtration of a simplicial complex is + * a function \f$f:\mathbf{K} \rightarrow \mathbb{R}\f$ satisfying \f$f(\tau)\leq f(\sigma)\f$ whenever + * \f$\tau \subseteq \sigma\f$. Ordering the simplices by increasing filtration values + * (breaking ties so as a simplex appears after its subsimplices of same filtration value) + * provides an indexing scheme. + * + + <DT>Implementations:</DT> + There are two implementation of complexes. The first on is the Simplex_tree data structure. + The simplex tree is an efficient and flexible + data structure for representing general (filtered) simplicial complexes. The data structure + is described in \cite boissonnatmariasimplextreealgorithmica + + The second one is the Hasse_complex. The Hasse complex is a data structure representing + explicitly all co-dimension 1 incidence relations in a complex. It is consequently faster + when accessing the boundary of a simplex, but is less compact and harder to construct from + scratch. + + + * \author Clément Maria + * \version 1.0 + * \date 2014 + * \copyright GNU General Public License v3. + * @{ + */ +/** + * \brief Simplex Tree data structure for representing simplicial complexes. + * + * \details Every simplex \f$[v_0, \cdots ,v_d]\f$ admits a canonical orientation + * induced by the order relation on vertices \f$ v_0 < \cdots < v_d \f$. + * + * Details may be found in \cite boissonnatmariasimplextreealgorithmica. + * + * \implements FilteredComplex + * + */ +template < typename IndexingTag = linear_indexing_tag + , typename FiltrationValue = double + , typename SimplexKey = int //must be a signed integer type + , typename VertexHandle = int //must be a signed integer type, int convertible to it +// , bool ContiguousVertexHandles = true //true is Vertex_handles are exactly the set [0;n) + > +class Simplex_tree +{ + +public: + typedef IndexingTag Indexing_tag; +/** \brief Type for the value of the filtration function. + * + * Must be comparable with <. */ + typedef FiltrationValue Filtration_value; +/** \brief Key associated to each simplex. + * + * Must be a signed integer type. */ + typedef SimplexKey Simplex_key; +/** \brief Type for the vertex handle. + * + * Must be a signed integer type. It admits a total order <. */ + typedef VertexHandle Vertex_handle; + + + /* Type of node in the simplex tree. */ + typedef Simplex_tree_node_explicit_storage < Simplex_tree > Node; + /* Type of dictionary Vertex_handle -> Node for traversing the simplex tree. */ + typedef typename boost::container::flat_map< Vertex_handle + , Node > Dictionary; + + + friend class Simplex_tree_node_explicit_storage < Simplex_tree < FiltrationValue + , SimplexKey + , VertexHandle > >; + friend class Simplex_tree_siblings < Simplex_tree < FiltrationValue + , SimplexKey + , VertexHandle > + , Dictionary > ; + friend class Simplex_tree_simplex_vertex_iterator < Simplex_tree < FiltrationValue + , SimplexKey + , VertexHandle > >; + friend class Simplex_tree_boundary_simplex_iterator < Simplex_tree < FiltrationValue + , SimplexKey + , VertexHandle > >; + friend class Simplex_tree_complex_simplex_iterator < Simplex_tree < FiltrationValue + , SimplexKey + , VertexHandle > >; + friend class Simplex_tree_skeleton_simplex_iterator < Simplex_tree < FiltrationValue + , SimplexKey + , VertexHandle > >; + template < class T1, class T2 > friend class Persistent_cohomology; + + /* \brief Set of nodes sharing a same parent in the simplex tree. */ + typedef Simplex_tree_siblings < Simplex_tree + , Dictionary > Siblings; +public: +/** \brief Handle type to a simplex contained in the simplicial complex represented + * byt he simplex tree. */ + typedef typename Dictionary::iterator Simplex_handle; +private: + typedef typename Dictionary::iterator Dictionary_it; + typedef typename Dictionary_it::value_type Dit_value_t; + + struct return_first { + Vertex_handle operator()(const Dit_value_t& p_sh) const {return p_sh.first;} + }; +public: +/** \name Range and iterator types + * + * The naming convention is Container_content_(iterator/range). A Container_content_range is + * essentially an object on which the methods begin() and end() can be called. They both return + * an object of type Container_content_iterator, and allow the traversal of the range + * [ begin();end() ). + * @{ */ + + /** \brief Iterator over the vertices of the simplicial complex. + * + * 'value_type' is Vertex_handle. */ + typedef boost::transform_iterator < return_first, Dictionary_it > Complex_vertex_iterator ; + /** \brief Range over the vertices of the simplicial complex. */ + typedef boost::iterator_range < Complex_vertex_iterator > Complex_vertex_range ; + /** \brief Iterator over the vertices of a simplex. + * + * 'value_type' is Vertex_handle. */ + typedef Simplex_tree_simplex_vertex_iterator < Simplex_tree > Simplex_vertex_iterator ; + /** \brief Range over the vertices of a simplex. */ + typedef boost::iterator_range < Simplex_vertex_iterator > Simplex_vertex_range ; + /** \brief Iterator over the simplices of the boundary of a simplex. + * + * 'value_type' is Simplex_handle. */ + 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 simplicial complex. + * + * 'value_type' is Simplex_handle. */ + typedef Simplex_tree_complex_simplex_iterator < Simplex_tree > Complex_simplex_iterator ; + /** \brief Range over the simplices of the simplicial complex. */ + typedef boost::iterator_range < Complex_simplex_iterator > Complex_simplex_range ; + /** \brief Iterator over the simplices of the skeleton of the simplicial complex, for a given + * dimension. + * + * 'value_type' is Simplex_handle. */ + typedef Simplex_tree_skeleton_simplex_iterator < Simplex_tree > Skeleton_simplex_iterator ; + /** \brief Range over the simplices of the skeleton of the simplicial complex, for a given + * dimension. */ + typedef boost::iterator_range < Skeleton_simplex_iterator > Skeleton_simplex_range ; + /** \brief Iterator over the simplices of the simplicial complex, ordered by the filtration. + * + * 'value_type' is Simplex_handle. */ + typedef typename std::vector < Simplex_handle >::iterator Filtration_simplex_iterator; + /** \brief Range over the simplices of the simplicial complex, ordered by the filtration. */ + typedef boost::iterator_range < Filtration_simplex_iterator > Filtration_simplex_range ; + +/* @} */ //end name range and iterator types + + +/** \name Range and iterator methods + * @{ */ + +/** \brief Returns a range over the vertices of the simplicial complex. + * + * The order is increasing according to < on Vertex_handles.*/ +Complex_vertex_range complex_vertex_range() +{ return Complex_vertex_range( boost::make_transform_iterator(root_.members_.begin(),return_first()) + , boost::make_transform_iterator(root_.members_.end(),return_first())); } +/** \brief Returns a range over the simplices of the simplicial complex. + * + * In the Simplex_tree, the tree is traverse in a depth-first fashion. + * Consequently, simplices are ordered according to lexicographic order on the list of + * Vertex_handles of a simplex, read in increasing < order for Vertex_handles. */ + Complex_simplex_range complex_simplex_range() + { return Complex_simplex_range ( Complex_simplex_iterator(this), + Complex_simplex_iterator() ); } +/** \brief Returns a range over the simplices of the dim-skeleton of the simplicial complex. + * + * The \f$d\f$-skeleton of a simplicial complex \f$\mathbf{K}\f$ is the simplicial complex containing the + * simplices of \f$\mathbf{K}\f$ of dimension at most \f$d\f$. + * + * @param[in] dim The maximal dimension of the simplices in the skeleton. + * + * The simplices are ordered according to lexicographic order on the list of + * Vertex_handles of a simplex, read in increasing < order for Vertex_handles. */ + Skeleton_simplex_range skeleton_simplex_range(int dim) + { return Skeleton_simplex_range ( Skeleton_simplex_iterator(this,dim) + , Skeleton_simplex_iterator() ); } +/** \brief Returns a range over the simplices of the simplicial complex, + * in the order of the filtration. + * + * The filtration is a monotonic function \f$ f: \mathbf{K} \rightarrow \mathbb{R} \f$, i.e. if two simplices + * \f$\tau\f$ and \f$\sigma\f$ satisfy \f$\tau \subseteq \sigma\f$ then + * \f$f(\tau) \leq f(\sigma)\f$. + * + * The method returns simplices ordered according to increasing filtration values. Ties are + * resolved by considering inclusion relation (subsimplices appear before their cofaces). If two + * simplices have same filtration value but are not comparable w.r.t. inclusion, lexicographic + * order is used. + * + * The filtration must be valid. If the filtration has not been initialized yet, the + * method initializes it (i.e. order the simplices). */ + Filtration_simplex_range filtration_simplex_range(linear_indexing_tag) + { if(filtration_vect_.empty()) { initialize_filtration(); } + return Filtration_simplex_range ( filtration_vect_.begin() + , filtration_vect_.end()); } + + Filtration_simplex_range filtration_simplex_range() { + return filtration_simplex_range(Indexing_tag()); + } +/** \brief Returns a range over the vertices of a simplex. + * + * The order in which the vertices are visited is the decreasing order for < on Vertex_handles, + * 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) + { return Simplex_vertex_range ( Simplex_vertex_iterator(this,sh), + Simplex_vertex_iterator(this));} + +/** \brief Returns a range over the simplices of the boundary of a simplex. + * + * 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$0\f$ to \f$d\f$, + * where \f$\widehat{v_i}\f$ means that the vertex \f$v_i\f$ is omitted. + * + * We note that the alternate sum of the simplices given by the iterator + * gives \f$(-1)^{\text{dim} \sigma}\f$ the chains corresponding to the boundary + * of the simplex. + * + * @param[in] sh Simplex for which the boundary is computed. */ + Boundary_simplex_range boundary_simplex_range(Simplex_handle sh) + { return Boundary_simplex_range ( Boundary_simplex_iterator(this,sh), + Boundary_simplex_iterator(this) ); } + +/** @} */ //end range and iterator methods + + +/** \name Constructor/Destructor + * @{ */ + +/** \brief Constructs an empty simplex tree. */ + Simplex_tree () + : null_vertex_(-1) + , threshold_(0) + , num_simplices_(0) + , root_(NULL,null_vertex_) + , filtration_vect_() + , dimension_(-1) {} + +/** \brief Destructor; deallocates the whole tree structure. */ +~Simplex_tree() +{ + for(auto sh = root_.members().begin(); sh != root_.members().end(); ++sh) + { + if(has_children(sh)) { rec_delete(sh->second.children()); } + } +} +/** @} */ // end constructor/destructor +private: +/** Recursive deletion. */ +void rec_delete(Siblings * sib) +{ + for(auto sh = sib->members().begin(); sh != sib->members().end(); ++sh) + { if(has_children(sh)) { rec_delete(sh->second.children()); } } + delete sib; +} + + +public: +/** \brief Returns the key associated to a simplex. + * + * The filtration must be initialized. */ +Simplex_key key ( Simplex_handle sh ) { return sh->second.key(); } +/** \brief Returns the simplex associated to a key. + * + * The filtration must be initialized. */ +Simplex_handle simplex ( Simplex_key key ) { return filtration_vect_[key]; } +/** \brief Returns the filtration value of a simplex. + * + * Called on the null_simplex, returns INFINITY. */ +Filtration_value filtration(Simplex_handle sh) +{ + if(sh != null_simplex()) { return sh->second.filtration(); } + else { return INFINITY; }//filtration(); } +} +/** \brief Returns an upper bound of the filtration values of the simplices. */ +Filtration_value filtration() +{ return threshold_; } +/** \brief Returns a Simplex_handle different from all Simplex_handles + * associated to the simplices in the simplicial complex. + * + * One can call filtration(null_simplex()). */ +Simplex_handle null_simplex() { return Dictionary_it(NULL); } +/** \brief Returns a key different for all keys associated to the + * simplices of the simplicial complex. */ +Simplex_key null_key() { return -1; } +/** \brief Returns a Vertex_handle different from all Vertex_handles associated + * to the vertices of the simplicial complex. */ +Vertex_handle null_vertex() { return null_vertex_; } +/** \brief Returns the number of vertices in the complex. */ +size_t num_vertices() { return root_.members_.size(); } +/** \brief Returns the number of simplices in the complex. + * + * Does not count the empty simplex. */ +size_t num_simplices() { return num_simplices_; } + +/** \brief Returns the dimension of a simplex. + * + * Must be different from null_simplex().*/ +int dimension(Simplex_handle sh) +{ Siblings * curr_sib = self_siblings(sh); + int dim = 0; + while(curr_sib != NULL) { ++dim; curr_sib = curr_sib->oncles(); } + return dim-1; +} +/** \brief Returns an upper bound on the dimension of the simplicial complex. */ +int dimension() +{ return dimension_; } + +/** \brief Returns true iff the node in the simplex tree pointed by + * sh has children.*/ +bool has_children(Simplex_handle sh) +{ return (sh->second.children()->parent() == sh->first); } + +/** \brief Given a range of Vertex_handles, returns the Simplex_handle + * of the simplex in the simplicial complex containing the corresponding + * vertices. Return null_simplex() if the simplex is not in the complex. + * + * The type RandomAccessVertexRange must be a range for which .begin() and + * .end() return random access iterators, with <CODE>value_type</CODE> + * <CODE>Vertex_handle</CODE>. + */ +template <class RandomAccessVertexRange > +Simplex_handle find(RandomAccessVertexRange & s) +{ + if(s.begin() == s.end()) std::cerr << "Empty simplex \n"; + + sort(s.begin(),s.end()); + + Siblings * tmp_sib = &root_; + Dictionary_it tmp_dit; + Vertex_handle last = s[s.size()-1]; + for(auto v : s) { + tmp_dit = tmp_sib->members_.find(v); + if(tmp_dit == tmp_sib->members_.end()) { return null_simplex(); } + if( !has_children(tmp_dit) && v != last) { return null_simplex(); } + tmp_sib = tmp_dit->second.children(); + } + return tmp_dit; +} + +/** \brief Returns the Simplex_handle corresponding to the 0-simplex + * representing the vertex with Vertex_handle v. */ +Simplex_handle find_vertex(Vertex_handle v) +{ return root_.members_.begin()+v; } +//{ return root_.members_.find(v); } + + +/** \brief Insert a simplex, represented by a range of Vertex_handles, in the simplicial complex. + * + * @param[in] simplex range of Vertex_handles, representing the vertices of the new simplex + * @param[in] filtration the filtration value assigned to the new simplex. + * The return type is a pair. 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 + * to the new simplex. + * If the insertion fails (the simplex is already there), the bool is set to false. If the insertion + * fails and the simplex already in the complex has a filtration value strictly bigger than 'filtration', + * we assign this simplex with the new value 'filtration', and set the Simplex_handle filed of the + * output pair to the Simplex_handle of the simplex. Otherwise, we set the Simplex_handle part to + * null_simplex. + * + * All subsimplices do not necessary need to be already in the simplex tree to proceed to an + * insertion. However, the property of being a simplicial complex will be violated. This allows + * us to insert a stream of simplices contained in a simplicial complex without considering any + * order on them. + * + * The filtration value + * assigned to the new simplex must preserve the monotonicity of the filtration. + * + * The type RandomAccessVertexRange must be a range for which .begin() and + * .end() return random access iterators, with 'value_type' Vertex_handle. */ +template <class RandomAccessVertexRange > +std::pair< Simplex_handle, bool > +insert ( RandomAccessVertexRange & simplex + , Filtration_value filtration ) +{ + if(simplex.empty()) { return std::pair< Simplex_handle, bool >(null_simplex(),true); } + + sort(simplex.begin(),simplex.end()); //must be sorted in increasing order + + Siblings * curr_sib = &root_; + std::pair< Simplex_handle, bool > res_insert; + typename RandomAccessVertexRange::iterator vi; + for( vi = simplex.begin(); vi != simplex.end()-1; ++vi ) + { + res_insert = curr_sib->members_.emplace(*vi, Node(curr_sib,filtration)); + if(!(has_children(res_insert.first))) + { res_insert.first->second.assign_children( new Siblings(curr_sib, *vi)); } + curr_sib = res_insert.first->second.children(); + } + res_insert = curr_sib->members_.emplace(*vi, Node(curr_sib,filtration)); + if(!res_insert.second) //if already in the complex + { + if(res_insert.first->second.filtration() > filtration) //if filtration value modified + { + res_insert.first->second.assign_filtration(filtration); + return res_insert; + } + return std::pair< Simplex_handle, bool > (null_simplex(),false);// if filtration value unchanged + } + //otherwise the insertion has succeeded + return res_insert; +} + +/** \brief Assign a value 'key' to the key of the simplex + * represented by the Simplex_handle 'sh'. */ +void assign_key(Simplex_handle sh, Simplex_key key) +{ sh->second.assign_key(key);} + +public: +/** Returns the two Simplex_handle corresponding to the endpoints of + * and edge. sh must point to a 1-dimensional simplex. This is an + * optimized version of the boundary computation. */ +std::pair<Simplex_handle,Simplex_handle> endpoints(Simplex_handle sh) +{ return std::pair<Simplex_handle,Simplex_handle>(root_.members_.find(sh->first) + , root_.members_.find(self_siblings(sh)->parent()) ); } + +/** Returns the Siblings containing a simplex.*/ +Siblings * self_siblings(Simplex_handle sh) +{ if(sh->second.children()->parent() == sh->first) return sh->second.children()->oncles(); + else return sh->second.children(); } + +// void display_simplex(Simplex_handle sh) +// { +// std::cout << " " << "[" << filtration(sh) << "] "; +// for( auto vertex : simplex_vertex_range(sh) ) +// { std::cout << vertex << " "; } +// } + + // void print(Simplex_handle sh, std::ostream& os = std::cout) + // { for(auto v : simplex_vertex_range(sh)) {os << v << " ";} + // os << std::endl; } + +public: +/** Returns a pointer to the root nodes of the simplex tree. */ +Siblings * root() { return &root_; } + + + + +public: +/** Set an upper bound for the filtration values. */ +void set_filtration(Filtration_value fil) +{ threshold_ = fil; } +/** Set a number of simplices for the simplicial complex. */ +void set_num_simplices(size_t num_simplices) +{ num_simplices_ = num_simplices; } +/** Set a dimension for the simplicial complex. */ +void set_dimension(int dimension) +{ dimension_ = dimension; } + +public: +/** \brief Initializes the filtrations, i.e. sort the + * simplices according to their order in the filtration and initializes all Simplex_keys. + * + * 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). + * + * The use of a depth-first traversal of the simplex tree, provided by + * complex_simplex_range(), combined with + * a stable sort is meant to optimize the order of simplices with same + * filtration value. The heuristic consists in inserting the cofaces of a + * simplex as soon as possible. + * + * Will be automatically called when calling filtration_simplex_range() + * if the filtration has never been initialized yet. */ +void initialize_filtration() +{ + filtration_vect_.clear(); + filtration_vect_.reserve(num_simplices()); + for(auto cpx_it = complex_simplex_range().begin(); + cpx_it != complex_simplex_range().end(); ++cpx_it) { filtration_vect_.push_back(*cpx_it); } + +//the stable sort ensures the ordering heuristic + std::stable_sort(filtration_vect_.begin(),filtration_vect_.end(),is_before_in_filtration(this)); +} + +private: +/** \brief Returns true iff the list of vertices of sh1 + * is smaller than the list of vertices of sh2 w.r.t. + * lexicographic order on the lists read in reverse. + * + * It defines a StrictWeakOrdering on simplices. The Simplex_vertex_iterators + * must traverse the Vertex_handle in decreasing order. Reverse lexicographic order satisfy + * the property that a subsimplex of a simplex is always strictly smaller with this order. */ +bool reverse_lexicographic_order(Simplex_handle sh1, Simplex_handle sh2) +{ + Simplex_vertex_range rg1 = simplex_vertex_range(sh1); + Simplex_vertex_range rg2 = simplex_vertex_range(sh2); + Simplex_vertex_iterator it1 = rg1.begin(); + Simplex_vertex_iterator it2 = rg2.begin(); + while(it1 != rg1.end() && it2 != rg2.end()) + { + if(*it1 == *it2) {++it1; ++it2;} + else { return *it1 < *it2; } + } + return ( (it1 == rg1.end()) && (it2 != rg2.end()) ); +} +/** \brief StrictWeakOrdering, for the simplices, defined by the filtration. + * + * It corresponds to the partial order + * induced by the filtration values, with ties resolved using reverse lexicographic order. + * Reverse lexicographic order has the property to always consider the subsimplex of a simplex + * to be smaller. The filtration function must be monotonic. */ +struct is_before_in_filtration { + is_before_in_filtration(Simplex_tree * st) : st_(st) {} + + bool operator()( const Simplex_handle sh1, + const Simplex_handle sh2 ) const + { + if(st_->filtration(sh1) != st_->filtration(sh2)) + { return st_->filtration(sh1) < st_->filtration(sh2); } + + return st_->reverse_lexicographic_order(sh1,sh2); //is sh1 a proper subface of sh2 + } + + Simplex_tree * st_; +}; + +public: +/** \brief Inserts a 1-skeleton in an empty Simplex_tree. + * + * The Simplex_tree must contain no simplex when the method is + * called. + * + * Inserts all vertices and edges given by a OneSkeletonGraph. + * OneSkeletonGraph must be a model of boost::AdjacencyGraph, + * boost::EdgeListGraph and boost::PropertyGraph. + * + * The vertex filtration value is accessible through the property tag + * vertex_filtration_t. + * The edge filtration value is accessible through the property tag + * edge_filtration_t. + * + * boost::graph_traits<OneSkeletonGraph>::vertex_descriptor + * must be Vertex_handle. + * boost::graph_traits<OneSkeletonGraph>::directed_category + * must be undirected_tag. */ + template< class OneSkeletonGraph > + void insert_graph( OneSkeletonGraph & skel_graph ) + { + assert(num_simplices() == 0); //the simplex tree must be empty + + if(boost::num_vertices(skel_graph) == 0) { return; } + if(num_edges(skel_graph) == 0) { dimension_ = 0; } + else { dimension_ = 1; } + + num_simplices_ = boost::num_vertices(skel_graph) + boost::num_edges(skel_graph); + root_.members_.reserve(boost::num_vertices(skel_graph)); + + typename boost::graph_traits<OneSkeletonGraph>::vertex_iterator v_it, v_it_end; + for(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) ) ); + } + typename boost::graph_traits<OneSkeletonGraph>::edge_iterator e_it, e_it_end; + for(tie(e_it,e_it_end) = boost::edges(skel_graph); + e_it != e_it_end; ++e_it) + { + auto u = source( *e_it, skel_graph ); auto v = target( *e_it, skel_graph ); + if( u < v ) // count edges only once { std::swap(u,v); } // u < v + { + auto sh = find_vertex(u); + if(! has_children(sh) ) { sh->second.assign_children(new Siblings(&root_,sh->first)); } + + sh->second.children()->members().emplace ( v + , Node( sh->second.children() + , boost::get( edge_filtration_t() + , skel_graph + , *e_it) )); + } + } +} +/** \brief Expands the Simplex_tree containing only its one skeleton + * until dimension max_dim. + * + * The expanded simplicial complex until dimension \f$d\f$ + * attached to a graph \f$G\f$ is the maximal simplicial complex of + * dimension at most \f$d\f$ admitting the graph \f$G\f$ as \f$1\f$-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. */ +void expansion(int max_dim) +{ + dimension_ = max_dim; + for(Dictionary_it root_it = root_.members_.begin(); + root_it != root_.members_.end(); ++root_it) + { + if(has_children(root_it)) + { siblings_expansion(root_it->second.children(), max_dim-1); } + } + dimension_ = max_dim - dimension_; +} +private: +/** \brief Recursive expansion of the simplex tree.*/ +void siblings_expansion ( Siblings * siblings, //must contain elements + int k) +{ + if(dimension_ > k) { dimension_ = k; } + if(k == 0) return; + Dictionary_it next = siblings->members().begin(); ++next; + + static std::vector< std::pair<Vertex_handle , Node> > inter; // static, not thread-safe. + 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); + if(has_children(root_sh)) + { + intersection(inter, //output intersection + next, //begin + siblings->members().end(),//end + root_sh->second.children()->members().begin(), + root_sh->second.children()->members().end(), + s_h->second.filtration()); + if(inter.size() != 0) + { this->num_simplices_ += inter.size(); + Siblings * new_sib = new Siblings(siblings, //oncles + s_h->first, //parent + inter); //boost::container::ordered_unique_range_t + inter.clear(); + s_h->second.assign_children(new_sib); + siblings_expansion(new_sib,k-1); + } + else { s_h->second.assign_children(siblings); //ensure the children property + inter.clear();} + } + } +} +/** \brief Intersects Dictionary 1 [begin1;end1) with Dictionary 2 [begin2,end2) + * and assigns the maximal possible Filtration_value to the Nodes. */ +void intersection ( std::vector< std::pair< Vertex_handle, Node > > & intersection + , Dictionary_it begin1 + , Dictionary_it end1 + , Dictionary_it begin2 + , Dictionary_it end2 + , Filtration_value filtration ) +{ + if(begin1 == end1 || begin2 == end2) return;// 0; + while( true ) + { + if( begin1->first == begin2->first ) + { + intersection.push_back(std::pair< Vertex_handle, Node >( begin1->first, + Node( NULL + , maximum( begin1->second.filtration() + , begin2->second.filtration() + , filtration ) + ) + ) + ); + ++begin1; ++begin2; + if( begin1 == end1 || begin2 == end2 ) return; + } + else + { + if( begin1->first < begin2->first ) + { ++begin1; if(begin1 == end1) return; } + else { ++begin2; if(begin2 == end2) return;} + } + } +} +/** Maximum over 3 values.*/ +Filtration_value maximum( Filtration_value a, + Filtration_value b, + Filtration_value c ) +{ Filtration_value max = ( a < b ) ? b : a; + return ( ( max < c ) ? c : max ); } + +public: +/** \brief Write the hasse diagram of the simplicial complex in os. + * + * Each row in the file correspond to a simplex. A line is written: + * dim idx_1 ... idx_k fil where dim is the dimension of the simplex, + * idx_1 ... idx_k are the row index (starting from 0) of the simplices of the boundary + * of the simplex, and fil is its filtration value. */ +void print_hasse(std::ostream & os) +{ + os << num_simplices() << " " << std::endl; + for(auto sh : filtration_simplex_range(Indexing_tag())) + { + os << dimension(sh) << " "; + for(auto b_sh : boundary_simplex_range(sh)) { os << key(b_sh) << " "; } + os << filtration(sh) << " \n"; + } +} + +private: + Vertex_handle null_vertex_; +/** \brief Upper bound on the filtration values of the simplices.*/ + Filtration_value threshold_ ; +/** \brief Total number of simplices in the complex, without the empty simplex.*/ + size_t num_simplices_ ; +/** \brief Set of simplex tree Nodes representing the vertices.*/ + Siblings root_ ; +/** \brief Simplices ordered according to a filtration.*/ + std::vector< Simplex_handle > filtration_vect_ ; +/** \brief Upper bound on the dimension of the simplicial complex.*/ + int dimension_ ; +}; + +// Print a Simplex_tree in os. +template< typename T1, typename T2, typename T3 > +std::ostream& operator<< ( std::ostream & os + , Simplex_tree< T1, T2, T3 > & st ) +{ + for(auto sh : st.filtration_simplex_range()) + { + os << st.dimension(sh) << " "; + for(auto v : st.simplex_vertex_range(sh)) + { os << v << " ";} + os << st.filtration(sh) << "\n";// TODO: VR - why adding the key ?? not read ?? << " " << st.key(sh) << " \n"; + } + return os; +} +template< typename T1, typename T2, typename T3 > +std::istream& operator>> ( std::istream & is + , Simplex_tree< T1, T2, T3 > & st ) +{ + //assert(st.num_simplices() == 0); + + std::vector< typename Simplex_tree<T1,T2,T3>::Vertex_handle > simplex; + typename Simplex_tree<T1,T2,T3>::Filtration_value fil; + typename Simplex_tree<T1,T2,T3>::Filtration_value max_fil = 0; + int max_dim = -1; + size_t num_simplices = 0; + while(read_simplex( is, simplex, fil )) //read all simplices in the file as a list of vertices + { + ++num_simplices; + int dim = (int)simplex.size() - 1; // Warning : simplex_size needs to be casted in int - Can be 0 + if(max_dim < dim) { max_dim = dim;} + if(max_fil < fil) { max_fil = fil;} + st.insert(simplex,fil); //insert every simplex in the simplex tree + simplex.clear(); + } + st.set_num_simplices(num_simplices); + st.set_dimension(max_dim); + st.set_filtration(max_fil); + + return is; +} + +/** @} */ //end defgroup simplex_tree + +} // namespace GUDHI + +#endif // GUDHI_FLAG_SIMPLEX_TREE_H |