diff options
author | vrouvrea <vrouvrea@636b058d-ea47-450e-bf9e-a15bfbe3eedb> | 2018-02-02 15:45:06 +0000 |
---|---|---|
committer | vrouvrea <vrouvrea@636b058d-ea47-450e-bf9e-a15bfbe3eedb> | 2018-02-02 15:45:06 +0000 |
commit | 5ae31a73253f6b29f5cbf1e04af17acf62741d9e (patch) | |
tree | 512041b9bf0835f2520e5b6b2fde0f972998efc4 /src/Simplex_tree/include/gudhi/Simplex_tree.h | |
parent | 4a5332dfb88ec27157c82f77df87f588f7016736 (diff) | |
parent | 265484997185f3bf900744406206a2d64ca0a20d (diff) |
Merge last trunk modificat
git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/rips_complex_from_correlation_matrix@3212 636b058d-ea47-450e-bf9e-a15bfbe3eedb
Former-commit-id: ea6296547b67ad82f5983d4309493e9131da7dd0
Diffstat (limited to 'src/Simplex_tree/include/gudhi/Simplex_tree.h')
-rw-r--r-- | src/Simplex_tree/include/gudhi/Simplex_tree.h | 217 |
1 files changed, 120 insertions, 97 deletions
diff --git a/src/Simplex_tree/include/gudhi/Simplex_tree.h b/src/Simplex_tree/include/gudhi/Simplex_tree.h index 7da767cb..7456cb1f 100644 --- a/src/Simplex_tree/include/gudhi/Simplex_tree.h +++ b/src/Simplex_tree/include/gudhi/Simplex_tree.h @@ -49,6 +49,7 @@ #include <initializer_list> #include <algorithm> // for std::max #include <cstdint> // for std::uint32_t +#include <iterator> // for std::distance namespace Gudhi { @@ -106,8 +107,9 @@ class Simplex_tree { }; struct Key_simplex_base_dummy { Key_simplex_base_dummy() {} - void assign_key(Simplex_key) { } - Simplex_key key() const { assert(false); return -1; } + // Undefined so it will not link + void assign_key(Simplex_key); + Simplex_key key() const; }; typedef typename std::conditional<Options::store_key, Key_simplex_base_real, Key_simplex_base_dummy>::type Key_simplex_base; @@ -121,7 +123,7 @@ class Simplex_tree { }; struct Filtration_simplex_base_dummy { Filtration_simplex_base_dummy() {} - void assign_filtration(Filtration_value f) { assert(f == 0); } + void assign_filtration(Filtration_value GUDHI_CHECK_code(f)) { GUDHI_CHECK(f == 0, "filtration value specified for a complex that does not store them"); } Filtration_value filtration() const { return 0; } }; typedef typename std::conditional<Options::store_filtration, Filtration_simplex_base_real, @@ -391,13 +393,13 @@ class Simplex_tree { return sh->second.key(); } - /** \brief Returns the simplex associated to a key. + /** \brief Returns the simplex that has index idx in the filtration. * * The filtration must be initialized. * \pre SimplexTreeOptions::store_key */ - Simplex_handle simplex(Simplex_key key) const { - return filtration_vect_[key]; + Simplex_handle simplex(Simplex_key idx) const { + return filtration_vect_[idx]; } /** \brief Returns the filtration value of a simplex. @@ -500,6 +502,7 @@ class Simplex_tree { * sh has children.*/ template<class SimplexHandle> bool has_children(SimplexHandle sh) const { + // Here we rely on the root using null_vertex(), which cannot match any real vertex. return (sh->second.children()->parent() == sh->first); } @@ -529,18 +532,30 @@ class Simplex_tree { Simplex_handle find_simplex(const std::vector<Vertex_handle> & simplex) { Siblings * tmp_sib = &root_; Dictionary_it tmp_dit; - Vertex_handle last = simplex.back(); - for (auto v : simplex) { - tmp_dit = tmp_sib->members_.find(v); - if (tmp_dit == tmp_sib->members_.end()) { + auto vi = simplex.begin(); + if (Options::contiguous_vertices) { + // Equivalent to the first iteration of the normal loop + GUDHI_CHECK(contiguous_vertices(), "non-contiguous vertices"); + Vertex_handle v = *vi++; + if(v < 0 || v >= static_cast<Vertex_handle>(root_.members_.size())) return null_simplex(); - } - if (!has_children(tmp_dit) && v != last) { + tmp_dit = root_.members_.begin() + v; + if (vi == simplex.end()) + return tmp_dit; + if (!has_children(tmp_dit)) + return null_simplex(); + tmp_sib = tmp_dit->second.children(); + } + for (;;) { + tmp_dit = tmp_sib->members_.find(*vi++); + if (tmp_dit == tmp_sib->members_.end()) + return null_simplex(); + if (vi == simplex.end()) + return tmp_dit; + if (!has_children(tmp_dit)) return null_simplex(); - } tmp_sib = tmp_dit->second.children(); } - return tmp_dit; } /** \brief Returns the Simplex_handle corresponding to the 0-simplex @@ -584,12 +599,14 @@ class Simplex_tree { std::pair<Simplex_handle, bool> res_insert; auto vi = simplex.begin(); for (; vi != simplex.end() - 1; ++vi) { + GUDHI_CHECK(*vi != null_vertex(), "cannot use the dummy null_vertex() as a real vertex"); res_insert = curr_sib->members_.emplace(*vi, Node(curr_sib, filtration)); if (!(has_children(res_insert.first))) { res_insert.first->second.assign_children(new Siblings(curr_sib, *vi)); } curr_sib = res_insert.first->second.children(); } + GUDHI_CHECK(*vi != null_vertex(), "cannot use the dummy null_vertex() as a real vertex"); res_insert = curr_sib->members_.emplace(*vi, Node(curr_sib, filtration)); if (!res_insert.second) { // if already in the complex @@ -664,71 +681,67 @@ class Simplex_tree { */ template<class InputVertexRange = std::initializer_list<Vertex_handle>> std::pair<Simplex_handle, bool> insert_simplex_and_subfaces(const InputVertexRange& Nsimplex, - Filtration_value filtration = 0) { + 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); // ----->> + return { null_simplex(), true }; // ----->> // Copy before sorting - std::vector<Vertex_handle> copy(first, last); + thread_local std::vector<Vertex_handle> copy; + copy.clear(); + copy.insert(copy.end(), first, last); std::sort(std::begin(copy), std::end(copy)); + GUDHI_CHECK_code( + for (Vertex_handle v : copy) + GUDHI_CHECK(v != null_vertex(), "cannot use the dummy null_vertex() as a real vertex"); + ) - 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); + return insert_simplex_and_subfaces_sorted(copy, 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); - } - 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 (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\n"; - exit(-1); + /// Same as insert_simplex_and_subfaces but assumes that the range of vertices is sorted + template<class ForwardVertexRange = std::initializer_list<Vertex_handle>> + std::pair<Simplex_handle, bool> insert_simplex_and_subfaces_sorted(const ForwardVertexRange& Nsimplex, Filtration_value filt = 0) { + auto first = std::begin(Nsimplex); + auto last = std::end(Nsimplex); + if (first == last) + return { null_simplex(), true }; // FIXME: false would make more sense to me. + GUDHI_CHECK(std::is_sorted(first, last), "simplex vertices listed in unsorted order"); + // Update dimension if needed. We could wait to see if the insertion succeeds, but I doubt there is much to gain. + dimension_ = (std::max)(dimension_, static_cast<int>(std::distance(first, last)) - 1); + return rec_insert_simplex_and_subfaces_sorted(root(), first, last, filt); + } + // To insert {1,2,3,4}, we insert {2,3,4} twice, once at the root, and once below 1. + template<class ForwardVertexIterator> + std::pair<Simplex_handle, bool> rec_insert_simplex_and_subfaces_sorted(Siblings* sib, ForwardVertexIterator first, ForwardVertexIterator last, Filtration_value filt) { + // An alternative strategy would be: + // - try to find the complete simplex, if found (and low filtration) exit + // - insert all the vertices at once in sib + // - loop over those (new or not) simplices, with a recursive call(++first, last) + Vertex_handle vertex_one = *first; + auto&& dict = sib->members(); + auto insertion_result = dict.emplace(vertex_one, Node(sib, filt)); + Simplex_handle simplex_one = insertion_result.first; + bool one_is_new = insertion_result.second; + if (!one_is_new) { + if (filtration(simplex_one) > filt) { + assign_filtration(simplex_one, filt); + } else { + // FIXME: this interface makes no sense, and it doesn't seem to be tested. + insertion_result.first = null_simplex(); } - 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 { - std::cerr << "Simplex_tree::rec_insert_simplex_and_subfaces - Recursivity error\n"; - exit(-1); } - return insert_result; + if (++first == last) return insertion_result; + if (!has_children(simplex_one)) + // TODO: have special code here, we know we are building the whole subtree from scratch. + simplex_one->second.assign_children(new Siblings(sib, vertex_one)); + auto res = rec_insert_simplex_and_subfaces_sorted(simplex_one->second.children(), first, last, filt); + // No need to continue if the full simplex was already there with a low enough filtration value. + if (res.first != null_simplex()) rec_insert_simplex_and_subfaces_sorted(sib, first, last, filt); + return res; } public: @@ -761,8 +774,12 @@ class Simplex_tree { return &root_; } - /** Set a dimension for the simplicial complex. */ + /** \brief Set a dimension for the simplicial complex. + * \details This function must be used with caution because it disables dimension recomputation when required + * (this recomputation can be triggered by `remove_maximal_simplex()` or `prune_above_filtration()`). + */ void set_dimension(int dimension) { + dimension_to_be_lowered_ = false; dimension_ = dimension; } @@ -937,8 +954,9 @@ class Simplex_tree { * called. * * Inserts all vertices and edges given by a OneSkeletonGraph. - * OneSkeletonGraph must be a model of boost::AdjacencyGraph, - * boost::EdgeListGraph and boost::PropertyGraph. + * OneSkeletonGraph must be a model of + * <a href="http://www.boost.org/doc/libs/1_65_1/libs/graph/doc/EdgeListGraph.html">boost::EdgeListGraph</a> + * and <a href="http://www.boost.org/doc/libs/1_65_1/libs/graph/doc/PropertyGraph.html">boost::PropertyGraph</a>. * * The vertex filtration value is accessible through the property tag * vertex_filtration_t. @@ -948,7 +966,10 @@ class Simplex_tree { * boost::graph_traits<OneSkeletonGraph>::vertex_descriptor * must be Vertex_handle. * boost::graph_traits<OneSkeletonGraph>::directed_category - * must be undirected_tag. */ + * must be undirected_tag. + * + * If an edge appears with multiplicity, the function will arbitrarily pick + * one representative to read the filtration value. */ template<class OneSkeletonGraph> void insert_graph(const OneSkeletonGraph& skel_graph) { // the simplex tree must be empty @@ -979,18 +1000,22 @@ class Simplex_tree { ++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))); + if (u == v) throw "Self-loops are not simplicial"; + // We cannot skip edges with the wrong orientation and expect them to + // come a second time with the right orientation, that does not always + // happen in practice. emplace() should be a NOP when an element with the + // same key is already there, so seeing the same edge multiple times is + // ok. + // Should we actually forbid multiple edges? That would be consistent + // with rejecting self-loops. + if (v < u) std::swap(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))); } } @@ -1134,7 +1159,7 @@ class Simplex_tree { to_be_inserted=false; break; } - filt = std::max(filt, filtration(border_child)); + filt = (std::max)(filt, filtration(border_child)); } if (to_be_inserted) { intersection.emplace_back(next->first, Node(nullptr, filt)); @@ -1145,7 +1170,7 @@ class Simplex_tree { Siblings * new_sib = new Siblings(siblings, // oncles simplex->first, // parent boost::adaptors::reverse(intersection)); // boost::container::ordered_unique_range_t - std::vector<Simplex_handle> blocked_new_sib_list; + std::vector<Vertex_handle> blocked_new_sib_vertex_list; // As all intersections are inserted, we can call the blocker function on all new_sib members for (auto new_sib_member = new_sib->members().begin(); new_sib_member != new_sib->members().end(); @@ -1153,17 +1178,19 @@ class Simplex_tree { bool blocker_result = block_simplex(new_sib_member); // new_sib member has been blocked by the blocker function // add it to the list to be removed - do not perform it while looping on it - if (blocker_result) - blocked_new_sib_list.push_back(new_sib_member); + if (blocker_result) { + blocked_new_sib_vertex_list.push_back(new_sib_member->first); + } } - bool removed = false; - for (auto& blocked_new_sib_member : blocked_new_sib_list){ - removed = removed || remove_maximal_simplex(blocked_new_sib_member); - } - if (removed) { + if (blocked_new_sib_vertex_list.size() == new_sib->members().size()) { + // Specific case where all have to be deleted + delete new_sib; // ensure the children property simplex->second.assign_children(siblings); } else { + for (auto& blocked_new_sib_member : blocked_new_sib_vertex_list) { + new_sib->members().erase(blocked_new_sib_member); + } // ensure recursive call simplex->second.assign_children(new_sib); siblings_expansion_with_blockers(new_sib, max_dim, k - 1, block_simplex); @@ -1328,7 +1355,7 @@ class Simplex_tree { if (sh_dimension >= dimension_) // Stop browsing as soon as the dimension is reached, no need to go furter return false; - new_dimension = std::max(new_dimension, sh_dimension); + new_dimension = (std::max)(new_dimension, sh_dimension); } dimension_ = new_dimension; return true; @@ -1338,16 +1365,14 @@ class Simplex_tree { public: /** \brief Remove a maximal simplex. * @param[in] sh Simplex handle on the maximal simplex to remove. - * @return a boolean value that is an implementation detail, and that the user is supposed to ignore * \pre Please check the simplex has no coface before removing it. * \exception std::invalid_argument In debug mode, if sh has children. - * \post Be aware that removing is shifting data in a flat_map (`initialize_filtration()` to be done). + * \post Be aware that removing is shifting data in a flat_map (initialize_filtration to be done). * \post Note that the dimension of the simplicial complex may be lower after calling `remove_maximal_simplex()` * than it was before. However, `upper_bound_dimension()` will return the old value, which remains a valid upper * bound. If you care, you can call `dimension()` to recompute the exact dimension. - * \internal @return true if the leaf's branch has no other leaves (branch's children has been re-assigned), false otherwise. */ - bool remove_maximal_simplex(Simplex_handle sh) { + void remove_maximal_simplex(Simplex_handle sh) { // Guarantee the simplex has no children GUDHI_CHECK(!has_children(sh), std::invalid_argument("Simplex_tree::remove_maximal_simplex - argument has children")); @@ -1365,9 +1390,7 @@ class Simplex_tree { delete child; // dimension may need to be lowered dimension_to_be_lowered_ = true; - return true; } - return false; } private: |