summaryrefslogtreecommitdiff
path: root/src/Simplex_tree/include/gudhi/Simplex_tree.h
diff options
context:
space:
mode:
Diffstat (limited to 'src/Simplex_tree/include/gudhi/Simplex_tree.h')
-rw-r--r--src/Simplex_tree/include/gudhi/Simplex_tree.h729
1 files changed, 498 insertions, 231 deletions
diff --git a/src/Simplex_tree/include/gudhi/Simplex_tree.h b/src/Simplex_tree/include/gudhi/Simplex_tree.h
index 5d3753ca..ea61fa2e 100644
--- a/src/Simplex_tree/include/gudhi/Simplex_tree.h
+++ b/src/Simplex_tree/include/gudhi/Simplex_tree.h
@@ -20,14 +20,17 @@
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
-#ifndef SRC_SIMPLEX_TREE_INCLUDE_GUDHI_SIMPLEX_TREE_H_
-#define SRC_SIMPLEX_TREE_INCLUDE_GUDHI_SIMPLEX_TREE_H_
+#ifndef SIMPLEX_TREE_H_
+#define SIMPLEX_TREE_H_
#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>
+#include <gudhi/reader_utils.h>
+#include <gudhi/graph_simplicial_complex.h>
+
#include <boost/container/flat_map.hpp>
#include <boost/iterator/transform_iterator.hpp>
#include <boost/graph/adjacency_list.hpp>
@@ -35,10 +38,11 @@
#include <algorithm>
#include <utility>
#include <vector>
-#include <iostream>
+#include <functional> // for greater<>
+#include <initializer_list>
namespace Gudhi {
-
+
/** \defgroup simplex_tree Filtered Complexes
*
* A simplicial complex \f$\mathbf{K}\f$
@@ -73,6 +77,17 @@ namespace Gudhi {
* \copyright GNU General Public License v3.
* @{
*/
+
+/// Model of SimplexTreeOptions.
+struct Simplex_tree_options_full_featured {
+ typedef linear_indexing_tag Indexing_tag;
+ typedef int Vertex_handle;
+ typedef double Filtration_value;
+ typedef int Simplex_key;
+ static const bool store_key = true;
+ static const bool store_filtration = true;
+};
+
/**
* \brief Simplex Tree data structure for representing simplicial complexes.
*
@@ -84,46 +99,69 @@ namespace Gudhi {
* \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)
->
+
+template<typename SimplexTreeOptions = Simplex_tree_options_full_featured>
class Simplex_tree {
public:
- typedef IndexingTag Indexing_tag;
+ typedef SimplexTreeOptions Options;
+ typedef typename Options::Indexing_tag Indexing_tag;
/** \brief Type for the value of the filtration function.
*
* Must be comparable with <. */
- typedef FiltrationValue Filtration_value;
+ typedef typename Options::Filtration_value Filtration_value;
/** \brief Key associated to each simplex.
*
* Must be a signed integer type. */
- typedef SimplexKey Simplex_key;
+ typedef typename Options::Simplex_key Simplex_key;
/** \brief Type for the vertex handle.
*
* Must be a signed integer type. It admits a total order <. */
- typedef VertexHandle Vertex_handle;
+ typedef typename Options::Vertex_handle 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. */
+ // Note: this wastes space when Vertex_handle is 32 bits and Node is aligned on 64 bits. It would be better to use a
+ // flat_set (with our own comparator) where we can control the layout of the struct (put Vertex_handle and
+ // Simplex_key next to each other).
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> >;
-
/* \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; }
+ Simplex_key key() const { return key_; }
+ private:
+ Simplex_key key_;
+ };
+ struct Key_simplex_base_dummy {
+ Key_simplex_base_dummy() {}
+ void assign_key(Simplex_key) { }
+ Simplex_key key() const { assert(false); return -1; }
+ };
+ typedef typename std::conditional<Options::store_key, Key_simplex_base_real, Key_simplex_base_dummy>::type 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_; }
+ private:
+ Filtration_value filt_;
+ };
+ struct Filtration_simplex_base_dummy {
+ Filtration_simplex_base_dummy() {}
+ void assign_filtration(Filtration_value f) { assert(f == 0); }
+ Filtration_value filtration() const { return 0; }
+ };
+ typedef typename std::conditional<Options::store_filtration, Filtration_simplex_base_real,
+ Filtration_simplex_base_dummy>::type Filtration_simplex_base;
+
public:
/** \brief Handle type to a simplex contained in the simplicial complex represented
- * byt he simplex tree. */
+ * by the simplex tree. */
typedef typename Dictionary::iterator Simplex_handle;
private:
@@ -157,6 +195,8 @@ class Simplex_tree {
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 Range over the cofaces of a simplex. */
+ typedef std::vector<Simplex_handle> Cofaces_simplex_range;
/** \brief Iterator over the simplices of the boundary of a simplex.
*
* 'value_type' is Simplex_handle. */
@@ -177,24 +217,23 @@ class Simplex_tree {
/** \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 Range over the simplices of the simplicial complex, ordered by the filtration. */
+ typedef std::vector<Simplex_handle> Filtration_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;
+ typedef typename Filtration_simplex_range::const_iterator Filtration_simplex_iterator;
/* @} */ // end name range and iterator types
/** \name Range and iterator methods
* @{ */
- /** \brief Returns a range over the vertices of the simplicial complex.
- *
+ /** \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()));
+ 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.
@@ -234,18 +273,15 @@ class Simplex_tree {
* 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(Indexing_tag) {
+ * 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. */
+ Filtration_simplex_range const& filtration_simplex_range(Indexing_tag = Indexing_tag()) {
if (filtration_vect_.empty()) {
initialize_filtration();
}
- return Filtration_simplex_range(filtration_vect_.begin(),
- filtration_vect_.end());
+ return filtration_vect_;
}
- 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,
@@ -253,6 +289,7 @@ class Simplex_tree {
* 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
return Simplex_vertex_range(Simplex_vertex_iterator(this, sh),
Simplex_vertex_iterator(this));
}
@@ -283,11 +320,48 @@ class Simplex_tree {
/** \brief Constructs an empty simplex tree. */
Simplex_tree()
: null_vertex_(-1),
- threshold_(0),
- num_simplices_(0),
- root_(NULL, null_vertex_),
- filtration_vect_(),
- dimension_(-1) {
+ threshold_(0),
+ root_(nullptr, null_vertex_),
+ filtration_vect_(),
+ dimension_(-1) { }
+
+ /** \brief User-defined copy constructor reproduces the whole tree structure. */
+ Simplex_tree(const Simplex_tree& simplex_source)
+ : null_vertex_(simplex_source.null_vertex_),
+ threshold_(simplex_source.threshold_),
+ filtration_vect_(),
+ dimension_(simplex_source.dimension_) {
+ auto root_source = simplex_source.root_;
+ auto memb_source = root_source.members();
+ root_ = Siblings(nullptr, null_vertex_, memb_source);
+ rec_copy(&root_, &root_source);
+ }
+
+ /** \brief depth first search, inserts simplices when reaching a leaf. */
+ void rec_copy(Siblings *sib, Siblings *sib_source) {
+ for (auto sh = sib->members().begin(), sh_source = sib_source->members().begin();
+ sh != sib->members().end(); ++sh, ++sh_source) {
+ if (has_children(sh_source)) {
+ Siblings * newsib = new Siblings(sib, sh_source->first);
+ newsib->members_.reserve(sh_source->second.children()->members().size());
+ for (auto & child : sh_source->second.children()->members())
+ newsib->members_.emplace_hint(newsib->members_.end(), child.first, Node(sib, child.second.filtration()));
+ rec_copy(newsib, sh_source->second.children());
+ sh->second.assign_children(newsib);
+ }
+ }
+ }
+
+ /** \brief User-defined move constructor moves the whole tree structure. */
+ Simplex_tree(Simplex_tree && old)
+ : null_vertex_(std::move(old.null_vertex_)),
+ threshold_(std::move(old.threshold_)),
+ root_(std::move(old.root_)),
+ filtration_vect_(std::move(old.filtration_vect_)),
+ dimension_(std::move(old.dimension_)) {
+ old.dimension_ = -1;
+ old.threshold_ = 0;
+ old.root_ = Siblings(nullptr, null_vertex_);
}
/** \brief Destructor; deallocates the whole tree structure. */
@@ -300,7 +374,7 @@ class Simplex_tree {
}
/** @} */ // end constructor/destructor
private:
- /** Recursive deletion. */
+ // Recursive deletion
void rec_delete(Siblings * sib) {
for (auto sh = sib->members().begin(); sh != sib->members().end(); ++sh) {
if (has_children(sh)) {
@@ -311,72 +385,135 @@ class Simplex_tree {
}
public:
+ /** \brief Checks if two simplex trees are equal. */
+ bool operator==(Simplex_tree& st2) {
+ if ((null_vertex_ != st2.null_vertex_) ||
+ (threshold_ != st2.threshold_) ||
+ (dimension_ != st2.dimension_))
+ return false;
+ return rec_equal(&root_, &st2.root_);
+ }
+
+ /** \brief Checks if two simplex trees are different. */
+ bool operator!=(Simplex_tree& st2) {
+ return (!(*this == st2));
+ }
+
+ private:
+ /** rec_equal: Checks recursively whether or not two simplex trees are equal, using depth first search. */
+ bool rec_equal(Siblings* s1, Siblings* s2) {
+ if (s1->members().size() != s2->members().size())
+ return false;
+ for (auto sh1 = s1->members().begin(), sh2 = s2->members().begin();
+ (sh1 != s1->members().end() && sh2 != s2->members().end()); ++sh1, ++sh2) {
+ if (sh1->first != sh2->first || sh1->second.filtration() != sh2->second.filtration())
+ return false;
+ if (has_children(sh1) != has_children(sh2))
+ return false;
+ // Recursivity on children only if both have children
+ else if (has_children(sh1))
+ if (!rec_equal(sh1->second.children(), sh2->second.children()))
+ return false;
+ }
+ return true;
+ }
+
+ public:
/** \brief Returns the key associated to a simplex.
*
- * The filtration must be initialized. */
- Simplex_key key(Simplex_handle sh) {
+ * The filtration must be initialized.
+ * \pre SimplexTreeOptions::store_key
+ */
+ static 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) {
+ * The filtration must be initialized.
+ * \pre SimplexTreeOptions::store_key
+ */
+ Simplex_handle simplex(Simplex_key key) const {
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) const {
+ * Called on the null_simplex, returns INFINITY.
+ * If SimplexTreeOptions::store_filtration is false, returns 0.
+ */
+ static 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() const {
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() const {
- return Dictionary_it(NULL);
+ static Simplex_handle null_simplex() {
+ return Dictionary_it(nullptr);
}
+
/** \brief Returns a key different for all keys associated to the
* simplices of the simplicial complex. */
- Simplex_key null_key() const {
+ static 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() const {
return null_vertex_;
}
+
/** \brief Returns the number of vertices in the complex. */
size_t num_vertices() const {
return root_.members_.size();
}
- /** \brief Returns the number of simplices in the complex.
- *
- * Does not count the empty simplex. */
- unsigned int num_simplices() const {
- return num_simplices_;
+
+ public:
+ /** \brief returns the number of simplices in the simplex_tree. */
+ size_t num_simplices() {
+ return num_simplices(&root_);
}
+ private:
+ /** \brief returns the number of simplices in the simplex_tree. */
+ size_t num_simplices(Siblings * sib) {
+ auto sib_begin = sib->members().begin();
+ auto sib_end = sib->members().end();
+ size_t simplices_number = sib_end - sib_begin;
+ for (auto sh = sib_begin; sh != sib_end; ++sh) {
+ if (has_children(sh)) {
+ simplices_number += num_simplices(sh->second.children());
+ }
+ }
+ return simplices_number;
+ }
+
+ public:
/** \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) {
+ while (curr_sib != nullptr) {
++dim;
curr_sib = curr_sib->oncles();
}
return dim - 1;
}
+
/** \brief Returns an upper bound on the dimension of the simplicial complex. */
int dimension() const {
return dimension_;
@@ -388,25 +525,34 @@ class Simplex_tree {
return (sh->second.children()->parent() == sh->first);
}
- /** \brief Given a range of Vertex_handles, returns the Simplex_handle
+ /** \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>.
+ * The type InputVertexRange must be a range of <CODE>Vertex_handle</CODE>
+ * on which we can call std::begin() function
*/
- template<class RandomAccessVertexRange>
- Simplex_handle find(RandomAccessVertexRange & s) {
- if (s.begin() == s.end())
- std::cerr << "Empty simplex \n";
-
- sort(s.begin(), s.end());
+ template<class InputVertexRange=std::initializer_list<Vertex_handle>>
+ Simplex_handle find(const InputVertexRange & s) {
+ auto first = std::begin(s);
+ auto last = std::end(s);
+
+ if (first == last)
+ return null_simplex(); // ----->>
+
+ // Copy before sorting
+ std::vector<Vertex_handle> copy(first, last);
+ std::sort(std::begin(copy), std::end(copy));
+ return find_simplex(copy);
+ }
+ private:
+ /** Find function, with a sorted range of vertices. */
+ Simplex_handle find_simplex(const std::vector<Vertex_handle> & simplex) {
Siblings * tmp_sib = &root_;
Dictionary_it tmp_dit;
- Vertex_handle last = s[s.size() - 1];
- for (auto v : s) {
+ Vertex_handle last = simplex.back();
+ for (auto v : simplex) {
tmp_dit = tmp_sib->members_.find(v);
if (tmp_dit == tmp_sib->members_.end()) {
return null_simplex();
@@ -424,8 +570,39 @@ class Simplex_tree {
Simplex_handle find_vertex(Vertex_handle v) {
return root_.members_.begin() + v;
}
-//{ return root_.members_.find(v); }
+ //{ return root_.members_.find(v); }
+
+ private:
+ /** \brief Inserts a simplex represented by a vector of vertex.
+ \warning the vector must be sorted by increasing vertex handle order */
+ std::pair<Simplex_handle, bool> insert_vertex_vector(const std::vector<Vertex_handle>& 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) {
+ 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;
+ }
+ // if filtration value unchanged
+ return std::pair<Simplex_handle, bool>(null_simplex(), false);
+ }
+ // otherwise the insertion has succeeded
+ return res_insert;
+ }
+ public:
/** \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
@@ -435,7 +612,7 @@ class Simplex_tree {
* 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
+ * we assign this simplex with the new value 'filtration', and set the Simplex_handle field of the
* output pair to the Simplex_handle of the simplex. Otherwise, we set the Simplex_handle part to
* null_simplex.
*
@@ -447,76 +624,107 @@ class Simplex_tree {
* 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_simplex(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;
+ * The type InputVertexRange must be a range for which .begin() and
+ * .end() return input iterators, with 'value_type' Vertex_handle. */
+ template<class InputVertexRange=std::initializer_list<Vertex_handle>>
+ std::pair<Simplex_handle, bool> insert_simplex(const InputVertexRange & simplex,
+ Filtration_value filtration = 0) {
+ auto first = std::begin(simplex);
+ auto last = std::end(simplex);
+
+ if (first == last)
+ return std::pair<Simplex_handle, bool>(null_simplex(), true); // ----->>
+
+ // Copy before sorting
+ std::vector<Vertex_handle> copy(first, last);
+ std::sort(std::begin(copy), std::end(copy));
+ return insert_vertex_vector(copy, filtration);
}
-
- /** \brief Insert a N-simplex and all his subfaces, from a N-simplex represented by a range of
+ /** \brief Insert a N-simplex and all his subfaces, from a N-simplex represented by a range of
* Vertex_handles, in the simplicial complex.
*
* @param[in] Nsimplex range of Vertex_handles, representing the vertices of the new N-simplex
* @param[in] filtration the filtration value assigned to the new N-simplex.
- */
- template<class RandomAccessVertexRange>
- void insert_simplex_and_subfaces(RandomAccessVertexRange& Nsimplex,
- Filtration_value filtration = 0.0) {
- if (Nsimplex.size() > 1) {
- for (unsigned int NIndex = 0; NIndex < Nsimplex.size(); NIndex++) {
- // insert N (N-1)-Simplex
- RandomAccessVertexRange NsimplexMinusOne;
- for (unsigned int NListIter = 0; NListIter < Nsimplex.size() - 1; NListIter++) {
- // (N-1)-Simplex creation
- NsimplexMinusOne.push_back( Nsimplex[(NIndex + NListIter) % Nsimplex.size()]);
- }
- // (N-1)-Simplex recursive call
- insert_simplex_and_subfaces(NsimplexMinusOne, filtration);
+ * 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 field of the
+ * output pair to the Simplex_handle of the simplex. Otherwise, we set the Simplex_handle part to
+ * null_simplex.
+ */
+ template<class InputVertexRange=std::initializer_list<Vertex_handle>>
+ std::pair<Simplex_handle, bool> insert_simplex_and_subfaces(const InputVertexRange& Nsimplex,
+ Filtration_value filtration = 0) {
+ auto first = std::begin(Nsimplex);
+ auto last = std::end(Nsimplex);
+
+ if (first == last)
+ return std::pair<Simplex_handle, bool>(null_simplex(), true); // ----->>
+
+ // Copy before sorting
+ std::vector<Vertex_handle> copy(first, last);
+ std::sort(std::begin(copy), std::end(copy));
+
+ std::vector<std::vector<Vertex_handle>> to_be_inserted;
+ std::vector<std::vector<Vertex_handle>> to_be_propagated;
+ return rec_insert_simplex_and_subfaces(copy, to_be_inserted, to_be_propagated, filtration);
+ }
+
+ private:
+ std::pair<Simplex_handle, bool> rec_insert_simplex_and_subfaces(std::vector<Vertex_handle>& the_simplex,
+ std::vector<std::vector<Vertex_handle>>& to_be_inserted,
+ std::vector<std::vector<Vertex_handle>>& to_be_propagated,
+ Filtration_value filtration = 0.0) {
+ std::pair<Simplex_handle, bool> insert_result;
+ if (the_simplex.size() > 1) {
+ // Get and remove last vertex
+ Vertex_handle last_vertex = the_simplex.back();
+ the_simplex.pop_back();
+ // Recursive call after last vertex removal
+ insert_result = rec_insert_simplex_and_subfaces(the_simplex, to_be_inserted, to_be_propagated, filtration);
+
+ // Concatenation of to_be_inserted and to_be_propagated
+ to_be_inserted.insert(to_be_inserted.begin(), to_be_propagated.begin(), to_be_propagated.end());
+ to_be_propagated = to_be_inserted;
+
+ // to_be_inserted treatment
+ for (auto& simplex_tbi : to_be_inserted) {
+ simplex_tbi.push_back(last_vertex);
}
- // N-Simplex insert
- std::pair<Simplex_handle, bool> returned = insert_simplex(Nsimplex, filtration);
- if (returned.second == true) {
- num_simplices_++;
+ std::vector<Vertex_handle> last_simplex(1, last_vertex);
+ to_be_inserted.insert(to_be_inserted.begin(), last_simplex);
+ // i.e. (0,1,2) =>
+ // [to_be_inserted | to_be_propagated] = [(1) (0,1) | (0)]
+ // [to_be_inserted | to_be_propagated] = [(2) (0,2) (1,2) (0,1,2) | (0) (1) (0,1)]
+ // N.B. : it is important the last inserted to be the highest in dimension
+ // in order to return the "last" insert_simplex result
+
+ // insert all to_be_inserted
+ for (auto& simplex_tbi : to_be_inserted) {
+ insert_result = insert_vertex_vector(simplex_tbi, filtration);
}
- } else if (Nsimplex.size() == 1) {
- // 1-Simplex insert - End of recursivity
- std::pair<Simplex_handle, bool> returned = insert_simplex(Nsimplex, filtration);
- if (returned.second == true) {
- num_simplices_++;
+ } else if (the_simplex.size() == 1) {
+ // When reaching the end of recursivity, vector of simplices shall be empty and filled on back recursive
+ if ((to_be_inserted.size() != 0) || (to_be_propagated.size() != 0)) {
+ std::cerr << "Simplex_tree::rec_insert_simplex_and_subfaces - Error vector not empty" << std::endl;
+ exit(-1);
}
+ std::vector<Vertex_handle> first_simplex(1, the_simplex.back());
+ // i.e. (0,1,2) => [to_be_inserted | to_be_propagated] = [(0) | ]
+ to_be_inserted.push_back(first_simplex);
+
+ insert_result = insert_vertex_vector(first_simplex, filtration);
} else {
- // Nothing to insert - empty vector
+ std::cerr << "Simplex_tree::rec_insert_simplex_and_subfaces - Recursivity error" << std::endl;
+ exit(-1);
}
+ return insert_result;
}
+ public:
/** \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) {
@@ -528,9 +736,8 @@ class Simplex_tree {
* 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()));
+ assert(dimension(sh) == 1);
+ return { find_vertex(sh->first), find_vertex(self_siblings(sh)->parent()) };
}
/** Returns the Siblings containing a simplex.*/
@@ -541,32 +748,17 @@ class Simplex_tree {
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(unsigned int num_simplices) {
- num_simplices_ = num_simplices;
- }
+
/** Set a dimension for the simplicial complex. */
void set_dimension(int dimension) {
dimension_ = dimension;
@@ -580,28 +772,114 @@ class Simplex_tree {
* 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
+ for (Simplex_handle sh : complex_simplex_range())
+ filtration_vect_.push_back(sh);
+
+ /* We use stable_sort here because with libstdc++ it is faster than sort.
+ * is_before_in_filtration is now a total order, but we used to call
+ * stable_sort for the following heuristic:
+ * 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.
+ */
std::stable_sort(filtration_vect_.begin(), filtration_vect_.end(),
is_before_in_filtration(this));
}
private:
+ /** Recursive search of cofaces
+ * This function uses DFS
+ *\param vertices contains a list of vertices, which represent the vertices of the simplex not found yet.
+ *\param curr_nbVertices represents the number of vertices of the simplex we reached by going through the tree.
+ *\param cofaces contains a list of Simplex_handle, representing all the cofaces asked.
+ *\param star true if we need the star of the simplex
+ *\param nbVertices number of vertices of the cofaces we search
+ * Prefix actions : When the bottom vertex matches with the current vertex in the tree, we remove the bottom vertex from vertices.
+ * Infix actions : Then we call or not the recursion.
+ * Postfix actions : Finally, we add back the removed vertex into vertices, and remove this vertex from curr_nbVertices so that we didn't change the parameters.
+ * If the vertices list is empty, we need to check if curr_nbVertices matches with the dimension of the cofaces asked.
+ */
+ void rec_coface(std::vector<Vertex_handle> &vertices, Siblings *curr_sib, int curr_nbVertices,
+ std::vector<Simplex_handle>& cofaces, bool star, int nbVertices) {
+ if (!(star || curr_nbVertices <= nbVertices)) // dimension of actual simplex <= nbVertices
+ return;
+ for (Simplex_handle simplex = curr_sib->members().begin(); simplex != curr_sib->members().end(); ++simplex) {
+ if (vertices.empty()) {
+ // 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
+ bool addCoface = (star || curr_nbVertices == nbVertices);
+ if (addCoface)
+ cofaces.push_back(simplex);
+ if ((!addCoface || star) && has_children(simplex)) // Rec call
+ rec_coface(vertices, simplex->second.children(), curr_nbVertices + 1, cofaces, star, nbVertices);
+ } else {
+ if (simplex->first == vertices.back()) {
+ // If curr_sib matches with the top vertex
+ bool equalDim = (star || curr_nbVertices == nbVertices); // dimension of actual simplex == nbVertices
+ bool addCoface = vertices.size() == 1 && equalDim;
+ if (addCoface)
+ cofaces.push_back(simplex);
+ if ((!addCoface || star) && has_children(simplex)) {
+ // Rec call
+ Vertex_handle tmp = vertices.back();
+ vertices.pop_back();
+ rec_coface(vertices, simplex->second.children(), curr_nbVertices + 1, cofaces, star, nbVertices);
+ vertices.push_back(tmp);
+ }
+ } else if (simplex->first > vertices.back()) {
+ return;
+ } else {
+ // (simplex->first < vertices.back()
+ if (has_children(simplex))
+ rec_coface(vertices, simplex->second.children(), curr_nbVertices + 1, cofaces, star, nbVertices);
+ }
+ }
+ }
+ }
+
+ public:
+ /** \brief Compute the star of a n simplex
+ * \param simplex represent the simplex of which we search the star
+ * \return Vector of Simplex_handle, empty vector if no cofaces found.
+ */
+
+ Cofaces_simplex_range star_simplex_range(const Simplex_handle simplex) {
+ return cofaces_simplex_range(simplex, 0);
+ }
+
+ /** \brief Compute the cofaces of a n simplex
+ * \param simplex represent the n-simplex of which we search the n+codimension cofaces
+ * \param codimension The function returns the n+codimension-cofaces of the n-simplex. If codimension = 0,
+ * return all cofaces (equivalent of star function)
+ * \return Vector of Simplex_handle, empty vector if no cofaces found.
+ */
+
+ Cofaces_simplex_range cofaces_simplex_range(const Simplex_handle simplex, int codimension) {
+ Cofaces_simplex_range cofaces;
+ // codimension must be positive or null integer
+ assert(codimension >= 0);
+ Simplex_vertex_range rg = simplex_vertex_range(simplex);
+ std::vector<Vertex_handle> copy(rg.begin(), rg.end());
+ if (codimension + static_cast<int>(copy.size()) > dimension_ + 1 ||
+ (codimension == 0 && static_cast<int>(copy.size()) > dimension_)) // n+codimension greater than dimension_
+ return cofaces;
+ // must be sorted in decreasing order
+ assert(std::is_sorted(copy.begin(), copy.end(), std::greater<Vertex_handle>()));
+ bool star = codimension == 0;
+ rec_coface(copy, &root_, 1, cofaces, star, codimension + static_cast<int>(copy.size()));
+ return cofaces;
+ }
+
+ 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.
@@ -624,6 +902,7 @@ class Simplex_tree {
}
return ((it1 == rg1.end()) && (it2 != rg2.end()));
}
+
/** \brief StrictWeakOrdering, for the simplices, defined by the filtration.
*
* It corresponds to the partial order
@@ -632,15 +911,15 @@ class Simplex_tree {
* to be smaller. The filtration function must be monotonic. */
struct is_before_in_filtration {
explicit is_before_in_filtration(Simplex_tree * st)
- : st_(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);
+ // Not using st_->filtration(sh1) because it uselessly tests for null_simplex.
+ if (sh1->second.filtration() != sh2->second.filtration()) {
+ return sh1->second.filtration() < sh2->second.filtration();
}
-
- return st_->reverse_lexicographic_order(sh1, sh2); // is sh1 a proper subface of sh2
+ // is sh1 a proper subface of sh2
+ return st_->reverse_lexicographic_order(sh1, sh2);
}
Simplex_tree * st_;
@@ -667,7 +946,8 @@ class Simplex_tree {
* must be undirected_tag. */
template<class OneSkeletonGraph>
void insert_graph(const OneSkeletonGraph& skel_graph) {
- assert(num_simplices() == 0); // the simplex tree must be empty
+ // the simplex tree must be empty
+ assert(num_simplices() == 0);
if (boost::num_vertices(skel_graph) == 0) {
return;
@@ -678,37 +958,37 @@ class Simplex_tree {
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 (std::tie(v_it, v_it_end) = boost::vertices(skel_graph); v_it != v_it_end;
- ++v_it) {
+ ++v_it) {
root_.members_.emplace_hint(
- root_.members_.end(), *v_it,
- Node(&root_, boost::get(vertex_filtration_t(), skel_graph, *v_it)));
+ 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 (std::tie(e_it, e_it_end) = boost::edges(skel_graph); e_it != e_it_end;
- ++e_it) {
+ ++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
+ 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)));
+ 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.
*
@@ -723,7 +1003,7 @@ class Simplex_tree {
void expansion(int max_dim) {
dimension_ = max_dim;
for (Dictionary_it root_it = root_.members_.begin();
- root_it != root_.members_.end(); ++root_it) {
+ root_it != root_.members_.end(); ++root_it) {
if (has_children(root_it)) {
siblings_expansion(root_it->second.children(), max_dim - 1);
}
@@ -734,7 +1014,7 @@ class Simplex_tree {
private:
/** \brief Recursive expansion of the simplex tree.*/
void siblings_expansion(Siblings * siblings, // must contain elements
- int k) {
+ int k) {
if (dimension_ > k) {
dimension_ = k;
}
@@ -745,68 +1025,55 @@ class Simplex_tree {
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) {
+ 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());
+ 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
+ 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
+ // ensure the children property
+ s_h->second.assign_children(siblings);
inter.clear();
}
}
}
}
+
/** \brief Intersects Dictionary 1 [begin1;end1) with Dictionary 2 [begin2,end2)
* and assigns the maximal possible Filtration_value to the Nodes. */
static 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) {
+ Filtration_value filtration_) {
if (begin1 == end1 || begin2 == end2)
return; // ----->>
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)
+ Filtration_value filt = (std::max)({begin1->second.filtration(), begin2->second.filtration(), filtration_});
+ intersection.emplace_back(begin1->first, Node(nullptr, filt));
+ if (++begin1 == end1 || ++begin2 == end2)
+ return; // ----->>
+ } else if (begin1->first < begin2->first) {
+ if (++begin1 == end1)
+ return;
+ } else /* begin1->first > begin2->first */ {
+ if (++begin2 == end2)
return; // ----->>
- } else {
- if (begin1->first < begin2->first) {
- ++begin1;
- if (begin1 == end1)
- return;
- } else {
- ++begin2;
- if (begin2 == end2)
- return; // ----->>
- }
}
}
}
- /** Maximum over 3 values.*/
- static 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.
@@ -817,7 +1084,7 @@ class Simplex_tree {
* 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())) {
+ for (auto sh : filtration_simplex_range()) {
os << dimension(sh) << " ";
for (auto b_sh : boundary_simplex_range(sh)) {
os << key(b_sh) << " ";
@@ -831,7 +1098,6 @@ class Simplex_tree {
/** \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.*/
- unsigned int num_simplices_;
/** \brief Set of simplex tree Nodes representing the vertices.*/
Siblings root_;
/** \brief Simplices ordered according to a filtration.*/
@@ -841,8 +1107,9 @@ class Simplex_tree {
};
// 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) {
+
+template<typename...T>
+std::ostream& operator<<(std::ostream & os, Simplex_tree<T...> & st) {
for (auto sh : st.filtration_simplex_range()) {
os << st.dimension(sh) << " ";
for (auto v : st.simplex_vertex_range(sh)) {
@@ -852,35 +1119,35 @@ std::ostream& operator<<(std::ostream & os, Simplex_tree<T1, T2, T3> & st) {
}
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;
+template<typename...T>
+std::istream& operator>>(std::istream & is, Simplex_tree<T...> & st) {
+ typedef Simplex_tree<T...> ST;
+ std::vector<typename ST::Vertex_handle> simplex;
+ typename ST::Filtration_value fil;
+ typename ST::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 = static_cast<int>(simplex.size() - 1); // Warning : simplex_size needs to be casted in int - Can be 0
+ while (read_simplex(is, simplex, fil)) {
+ // read all simplices in the file as a list of vertices
+ // Warning : simplex_size needs to be casted in int - Can be 0
+ int dim = static_cast<int> (simplex.size() - 1);
if (max_dim < dim) {
max_dim = dim;
}
if (max_fil < fil) {
max_fil = fil;
}
- st.insert_simplex(simplex, fil); // insert every simplex in the simplex tree
+ // insert every simplex in the simplex tree
+ st.insert_simplex(simplex, fil);
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 // SRC_SIMPLEX_TREE_INCLUDE_GUDHI_SIMPLEX_TREE_H_
+#endif // SIMPLEX_TREE_H_
generated by cgit v1.2.3 (git 2.39.1) at 2024-07-13 15:01:20 +0000