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.h179
1 files changed, 100 insertions, 79 deletions
diff --git a/src/Simplex_tree/include/gudhi/Simplex_tree.h b/src/Simplex_tree/include/gudhi/Simplex_tree.h
index cb6ab309..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,
@@ -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:
@@ -941,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.
@@ -952,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
@@ -983,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)));
}
}
@@ -1138,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));
@@ -1334,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;
generated by cgit v1.2.3 (git 2.39.1) at 2024-06-14 08:10:33 +0000