diff options
Diffstat (limited to 'src/Simplex_tree')
| -rw-r--r-- | src/Simplex_tree/include/gudhi/Simplex_tree.h | 296 | ||||
| -rw-r--r-- | src/Simplex_tree/test/simplex_tree_unit_test.cpp | 152 |
2 files changed, 290 insertions, 158 deletions
diff --git a/src/Simplex_tree/include/gudhi/Simplex_tree.h b/src/Simplex_tree/include/gudhi/Simplex_tree.h index d4a52113..6b016fc3 100644 --- a/src/Simplex_tree/include/gudhi/Simplex_tree.h +++ b/src/Simplex_tree/include/gudhi/Simplex_tree.h @@ -26,11 +26,13 @@ #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/reader_utils.h> #include <gudhi/Simplex_tree/indexing_tag.h> #include <boost/container/flat_map.hpp> #include <boost/iterator/transform_iterator.hpp> #include <boost/graph/adjacency_list.hpp> +#include <gudhi/graph_simplicial_complex.h> #include <algorithm> #include <utility> @@ -232,8 +234,7 @@ class Simplex_tree { Filtration_simplex_range filtration_simplex_range(Indexing_tag) { if (filtration_vect_.empty()) { initialize_filtration(); - } - return Filtration_simplex_range(filtration_vect_.begin(), + } return Filtration_simplex_range(filtration_vect_.begin(), filtration_vect_.end()); } @@ -284,6 +285,41 @@ class Simplex_tree { filtration_vect_(), dimension_(-1) { } + /** \brief Copy; copy the whole tree structure. */ + Simplex_tree(Simplex_tree& copy) : null_vertex_(copy.null_vertex_), + threshold_(copy.threshold_), + num_simplices_(copy.num_simplices_), + root_(NULL, -1, std::vector<std::pair<Vertex_handle, Node>> (copy.root_.members().begin(), copy.root_.members().end())), + filtration_vect_(copy.filtration_vect_), + dimension_(copy.dimension_) { + rec_copy(&root_, ©.root_); + } + + /** rec_copy: depth first search, inserts simplices when reaching a leaf.*/ + void rec_copy(Siblings *sib, Siblings *sib_copy) + { + for (auto sh = sib->members().begin(), sh_copy = sib_copy->members().begin(); sh != sib->members().end(); ++sh, ++sh_copy) { + if (has_children(sh_copy)) { + std::vector<std::pair<Vertex_handle, Node>> v(sh_copy->second.children()->members().begin(), sh_copy->second.children()->members().end()); + Siblings * newsib = new Siblings (sib, sh_copy->first); + for (auto it = v.begin(); it != v.end(); ++it) + newsib->members_.emplace(it->first, Node(sib, it->second.filtration())); + rec_copy(newsib, sh_copy->second.children()); + sh->second.assign_children(newsib); + } + } + } + + + + + /** \brief Move; moves the whole tree structure. */ + Simplex_tree(Simplex_tree&& old) : null_vertex_(std::move(old.null_vertex_)), threshold_(std::move(old.threshold_)), num_simplices_(std::move(old.num_simplices_)), root_(std::move(old.root_)), filtration_vect_(std::move(old.filtration_vect_)), dimension_(std::move(old.dimension_)) { + old.dimension_ = -1; + old.num_simplices_ = 0; + old.threshold_ = 0; + } + /** \brief Destructor; deallocates the whole tree structure. */ ~Simplex_tree() { for (auto sh = root_.members().begin(); sh != root_.members().end(); ++sh) { @@ -304,6 +340,72 @@ class Simplex_tree { delete sib; } + public: + /** \brief Prints the simplex_tree hierarchically. */ + void print_tree() + { + for (auto sh = root_.members().begin(); sh != root_.members().end(); ++sh) + { + std::cout << sh->first << " "; + if (has_children(sh)) + { + std::cout << "("; + rec_print(sh->second.children()); + std::cout << ")"; + } + std::cout << std::endl; + } + } + + + /** \brief Recursively prints the simplex_tree, using depth first search. */ + private: + void rec_print(Siblings * sib) + { + for (auto sh = sib->members().begin(); sh != sib->members().end(); ++sh) + { + std::cout << " " << sh->first << " "; + if (has_children(sh)) + { + std::cout << "("; + rec_print(sh->second.children()); + std::cout << ")"; + } + } + } + + public: + + /** \brief Checks whether or not two simplex trees are equal. */ + bool operator ==(Simplex_tree st2) + { + return (this->null_vertex_ == st2.null_vertex_ + && this->threshold_ == st2.threshold_ + && this->num_simplices_ == st2.num_simplices_ + && this->dimension_ == st2.dimension_ + && rec_equal(&root_, &st2.root_)); + } + + /** rec_equal: Checks recursively whether or not two simplex trees are equal, using depth first search. */ + private: + 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(); ++sh1, ++sh2) + { + if (sh1->first != sh2->first || sh1->second.filtration() != sh2->second.filtration()) + return false; + if (has_children(sh1)) + return rec_equal(sh1->second.children(), sh2->second.children()); + else if (has_children(sh2)) + return false; + else + return true; + } + return true; + } + public: /** \brief Returns the key associated to a simplex. * @@ -399,32 +501,36 @@ class Simplex_tree { return dimension_; } + /** \brief Returns true iff the node in the simplex tree pointed by * sh has children.*/ bool has_children(Simplex_handle sh) const { return (sh->second.children()->parent() == sh->first); } - public: /** \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()) // Empty simplex - return null_simplex(); - - sort(s.begin(), s.end()); +template<class InputVertexRange> + Simplex_handle find(const InputVertexRange & s) { + std::vector<Vertex_handle> copy(std::begin(s), std::end(s)); + 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(std::vector<Vertex_handle> & simplex) { + if (simplex.begin() == simplex.end()) + return null_simplex(); Siblings * tmp_sib = &root_; Dictionary_it tmp_dit; - Vertex_handle last = s[s.size() - 1]; - for (auto v : s) { + Vertex_handle last = simplex[simplex.size() - 1]; + for (auto v : simplex) { tmp_dit = tmp_sib->members_.find(v); if (tmp_dit == tmp_sib->members_.end()) { return null_simplex(); @@ -444,38 +550,14 @@ class Simplex_tree { } //{ 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_simplex(RandomAccessVertexRange & simplex, + private: + /** Recursively insert a simplex */ + template <typename RandomAccessVertexRange> + std::pair<Simplex_handle, bool> insert_simplex_rec(RandomAccessVertexRange & simplex, Filtration_value filtration) { if (simplex.empty()) { return std::pair<Simplex_handle, bool>(null_simplex(), true); } - // must be sorted in increasing order - sort(simplex.begin(), simplex.end()); - Siblings * curr_sib = &root_; std::pair<Simplex_handle, bool> res_insert; typename RandomAccessVertexRange::iterator vi; @@ -484,7 +566,7 @@ class Simplex_tree { if (!(has_children(res_insert.first))) { res_insert.first->second.assign_children(new Siblings(curr_sib, *vi)); } - curr_sib = res_insert.first->second.children(); + curr_sib = res_insert.first->second.children(); } res_insert = curr_sib->members_.emplace(*vi, Node(curr_sib, filtration)); if (!res_insert.second) { @@ -500,6 +582,37 @@ class Simplex_tree { // 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 + * @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_simplex(const RandomAccessVertexRange & simplex, + Filtration_value filtration) { + std::vector<Vertex_handle> copy(std::begin(simplex), std::end(simplex));; + std::sort(std::begin(copy), std::end(copy)); + return insert_simplex_rec(copy, filtration); +} /** \brief Insert a N-simplex and all his subfaces, from a N-simplex represented by a range of * Vertex_handles, in the simplicial complex. @@ -507,9 +620,21 @@ class Simplex_tree { * @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, + template<class InputVertexRange> + void insert_simplex_and_subfaces(const InputVertexRange& Nsimplex, Filtration_value filtration = 0.0) { + std::vector<Vertex_handle> copy(std::begin(Nsimplex), std::end(Nsimplex));; + std::sort(std::begin(copy), std::end(copy)); + insert_simplex_and_subfaces_rec(copy, filtration); + } + + private: + + /** Recursively insert simplex and all of its subfaces */ + template <typename RandomAccessVertexRange> + void insert_simplex_and_subfaces_rec(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 @@ -522,15 +647,17 @@ class Simplex_tree { insert_simplex_and_subfaces(NsimplexMinusOne, filtration); } // N-Simplex insert - std::pair<Simplex_handle, bool> returned = insert_simplex(Nsimplex, filtration); + std::pair<Simplex_handle, bool> returned = insert_simplex_rec(Nsimplex, filtration); } else if (Nsimplex.size() == 1) { // 1-Simplex insert - End of recursivity - std::pair<Simplex_handle, bool> returned = insert_simplex(Nsimplex, filtration); + std::pair<Simplex_handle, bool> returned = insert_simplex_rec(Nsimplex, filtration); } else { // Nothing to insert - empty vector } } + 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) { @@ -608,6 +735,81 @@ class Simplex_tree { is_before_in_filtration(this)); } +/** \brief Contracts two vertices : the contracted vertex is erased from the three, and the remaining vertex receives all the links of the contracted one. + \param remaining The value of the vertex within the other is contracted + \param deleted The value of the vertex to be contracted +*/ + void edge_contraction(Vertex_handle remaining, Vertex_handle deleted) + { + std::vector<Vertex_handle> accessRemaining, accessDeleted; + accessRemaining.push_back(remaining); + accessDeleted.push_back(deleted); + Simplex_handle shr = find(accessRemaining), shd = find(accessDeleted); + if (has_children(shd)) + { + Siblings * sibDeleted, * sibRemaining; // We need the siblings + sibDeleted = shd->second.children(); + sibRemaining = shr->second.children(); + rec_insert(sibDeleted, sibRemaining, shd->first); + } + rec_delete(&root_, shd->first); + } + + +/** \brief recursively insert the members of a Sibling into another, any time the replacedVertex appears into the target Sibling + \param sibInserted The sibling to be inserted + \param sibTarget The target sibling on which the other sibling is copied + \param replacedVertex The replacedVertex (Vertex_handle) needed to insert the sibling +*/ + void rec_insert(Siblings * sibInserted, Siblings * sibTarget, Vertex_handle replacedVertex) + { + for (auto sh = sibTarget->members().begin(); sh != sibTarget->members().end(); ++sh) + { + if (has_children(sh)) + rec_insert(sibInserted, sh->second.children(), replacedVertex); + if (sh->first == replacedVertex) + insert_members(sibTarget, sibInserted); + } + } + +/** \brief Copy a Sibling into another, which is possibly not empty + \param sibInserted The sibling to be copied + \param sibTarget The sibling on which we wan't to copy the other sibling +*/ + void insert_members(Siblings * sibTarget, Siblings * sibInserted) + { + for (auto sh = sibInserted->members().begin(); sh != sibInserted->members().end(); ++sh) + { + std::pair<Vertex_handle, Node> member(sh->first, Node(sibTarget, sh->second.filtration())); + if (has_children(sh)) + { + std::vector<std::pair<Vertex_handle, Node>> v(sh->second.children()->members().begin(), sh->second.children()->members().end()); + Siblings * newsib = new Siblings (sibTarget, sh->first); + for (auto it = v.begin(); it != v.end(); ++it) + newsib->members_.emplace(it->first, Node(sibTarget, it->second.filtration())); + insert_members(newsib, sh->second.children()); + member.second.assign_children(newsib); + + } + sibTarget->members_.emplace(member); + } + } + +/** \brief Erase every occurencies of a vertex in a Sibling + \param sib The sibling on which we wan't to erase the elements + \param deletedVertex The value of the members we wan't to erase +*/ + void rec_delete(Siblings * sib, Vertex_handle deletedVertex) + { + for (auto sh = sib->members().begin(); sh != sib->members().end(); ++sh) + { + if (has_children(sh)) + rec_delete(sh->second.children(), deletedVertex); + if (sh->first == deletedVertex) + sib->members_.erase(sh); + } + } + private: /** Recursive search of cofaces * This function uses DFS diff --git a/src/Simplex_tree/test/simplex_tree_unit_test.cpp b/src/Simplex_tree/test/simplex_tree_unit_test.cpp index 20403e2a..0c710651 100644 --- a/src/Simplex_tree/test/simplex_tree_unit_test.cpp +++ b/src/Simplex_tree/test/simplex_tree_unit_test.cpp @@ -7,7 +7,7 @@ #define BOOST_TEST_DYN_LINK #define BOOST_TEST_MODULE "simplex_tree" -#include <boost/test/unit_test.hpp> +#include <boost/test/included/unit_test.hpp> #include "gudhi/graph_simplicial_complex.h" #include "gudhi/reader_utils.h" @@ -125,8 +125,10 @@ void test_simplex_tree_contains(typeST& simplexTree, typeSimplex& simplex, int p std::cout << "test_simplex_tree_contains - filtration=" << simplexTree.filtration(*f_simplex) << "||" << simplex.second << std::endl; BOOST_CHECK(AreAlmostTheSame(simplexTree.filtration(*f_simplex), simplex.second)); - int simplexIndex = simplex.first.size() - 1; - for (auto vertex : simplexTree.simplex_vertex_range(*f_simplex)) { + int simplexIndex=simplex.first.size()-1; + std::sort(simplex.first.begin(), simplex.first.end()); // if the simplex wasn't sorted, the next test could fail + for( auto vertex : simplexTree.simplex_vertex_range(*f_simplex) ) + { std::cout << "test_simplex_tree_contains - vertex=" << vertex << "||" << simplex.first.at(simplexIndex) << std::endl; BOOST_CHECK(vertex == simplex.first.at(simplexIndex)); BOOST_CHECK(simplexIndex >= 0); @@ -521,113 +523,41 @@ BOOST_AUTO_TEST_CASE(NSimplexAndSubfaces_tree_insertion) { std::cout << std::endl; } - std::cout << "********************************************************************" << std::endl; - // TEST COFACE ALGORITHM - st.set_dimension(3); - std::cout << "COFACE ALGORITHM" << std::endl; - std::vector<Vertex_handle> v; - std::vector<Vertex_handle> simplex; - std::vector<typeST::Simplex_handle> result; - v.push_back(3); - std::cout << "First test : " << std::endl; - std::cout << "Star of (3):" << std::endl; - - simplex.push_back(3); - result.push_back(st.find(simplex)); - simplex.clear(); - - simplex.push_back(3); - simplex.push_back(0); - result.push_back(st.find(simplex)); - simplex.clear(); - - simplex.push_back(4); - simplex.push_back(3); - result.push_back(st.find(simplex)); - simplex.clear(); - - simplex.push_back(5); - simplex.push_back(4); - simplex.push_back(3); - result.push_back(st.find(simplex)); - simplex.clear(); - - simplex.push_back(5); - simplex.push_back(3); - result.push_back(st.find(simplex)); - simplex.clear(); - - test_cofaces(st, v, 0, result); - v.clear(); - result.clear(); - - v.push_back(1); - v.push_back(7); - std::cout << "Second test : " << std::endl; - std::cout << "Star of (1,7): " << std::endl; - - simplex.push_back(7); - simplex.push_back(1); - result.push_back(st.find(simplex)); - simplex.clear(); - - simplex.push_back(7); - simplex.push_back(6); - simplex.push_back(1); - simplex.push_back(0); - result.push_back(st.find(simplex)); - simplex.clear(); - - simplex.push_back(7); - simplex.push_back(1); - simplex.push_back(0); - result.push_back(st.find(simplex)); - simplex.clear(); - - simplex.push_back(7); - simplex.push_back(6); - simplex.push_back(1); - result.push_back(st.find(simplex)); - simplex.clear(); - - test_cofaces(st, v, 0, result); - result.clear(); - - std::cout << "Third test : " << std::endl; - std::cout << "2-dimension Cofaces of simplex(1,7) : " << std::endl; - - simplex.push_back(7); - simplex.push_back(1); - simplex.push_back(0); - result.push_back(st.find(simplex)); - simplex.clear(); - - simplex.push_back(7); - simplex.push_back(6); - simplex.push_back(1); - result.push_back(st.find(simplex)); - simplex.clear(); - - test_cofaces(st, v, 1, result); - result.clear(); - - std::cout << "Cofaces with a codimension too high (codimension + vetices > tree.dimension) :" << std::endl; - test_cofaces(st, v, 5, result); - // std::cout << "Cofaces with an empty codimension" << std::endl; - // test_cofaces(st, v, -1, result); - // std::cout << "Cofaces in an empty simplex tree" << std::endl; - // typeST empty_tree; - // test_cofaces(empty_tree, v, 1, result); - // std::cout << "Cofaces of an empty simplex" << std::endl; - // v.clear(); - // test_cofaces(st, v, 1, result); - - /* - // TEST Off read - std::cout << "********************************************************************" << std::endl; - typeST st2; - st2.tree_from_off("test.off"); - std::cout << st2; - */ + std::cout << "********************************************************************" << std::endl; + // TEST Copy constructor / Move + std::cout << "Printing st" << std::endl; + std::cout << &st << std::endl; + st.print_tree(); + typeST st_copy_1 = st; + BOOST_CHECK(st == st_copy_1); + typeST st_move = std::move(st); + std::cout << "Printing a copy of st" << std::endl; + std::cout << &st_copy_1 << std::endl; + st_copy_1.print_tree(); + BOOST_CHECK(st_move == st_copy_1); + std::cout << "Printing a move of st" << std::endl; + std::cout << &st_move << std::endl; + st_move.print_tree(); + typeST st_empty; + BOOST_CHECK(st == st_empty); + std::cout << "Printing st again" << std::endl; + std::cout << &st << std::endl; + st.print_tree(); + + std::cout << "********************************************************************" << std::endl; + // TEST Edge_contraction + typeST st_copy_2 = st_copy_1, st_copy_3 = st_copy_1, st_copy_4 = st_copy_1; + st_copy_1.edge_contraction(0, 3); + std::cout << "Printing a copy of st, with the edge (0, 3) contracted, 3 being contracted in 0" << std::endl; + st_copy_1.print_tree(); + st_copy_2.edge_contraction(1, 3); + std::cout << "Printing a copy of st, with the edge (1, 3) contracted, 3 being contracted in 1" << std::endl; + st_copy_2.print_tree(); + st_copy_3.edge_contraction(3, 4); + std::cout << "Printing a copy of st, with the edge (3, 4) contracted, 4 being contracted in 3" << std::endl; + st_copy_3.print_tree(); + st_copy_4.edge_contraction(1, 6); + std::cout << "Printing a copy of st, with the edge (1, 6) contracted, 6 being contracted in 1" << std::endl; + st_copy_4.print_tree(); } |