summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--src/Simplex_tree/include/gudhi/Simplex_tree.h71
-rw-r--r--src/Simplex_tree/test/simplex_tree_unit_test.cpp32
-rw-r--r--src/common/include/gudhi/distance_functions.h31
3 files changed, 103 insertions, 31 deletions
diff --git a/src/Simplex_tree/include/gudhi/Simplex_tree.h b/src/Simplex_tree/include/gudhi/Simplex_tree.h
index 580f85ed..7456cb1f 100644
--- a/src/Simplex_tree/include/gudhi/Simplex_tree.h
+++ b/src/Simplex_tree/include/gudhi/Simplex_tree.h
@@ -502,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);
}
@@ -531,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
@@ -586,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
@@ -678,6 +693,10 @@ class Simplex_tree {
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");
+ )
return insert_simplex_and_subfaces_sorted(copy, filtration);
}
@@ -935,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.
@@ -946,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
@@ -977,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)));
}
}
diff --git a/src/Simplex_tree/test/simplex_tree_unit_test.cpp b/src/Simplex_tree/test/simplex_tree_unit_test.cpp
index 580d610a..f63ea080 100644
--- a/src/Simplex_tree/test/simplex_tree_unit_test.cpp
+++ b/src/Simplex_tree/test/simplex_tree_unit_test.cpp
@@ -863,3 +863,35 @@ BOOST_AUTO_TEST_CASE(make_filtration_non_decreasing) {
BOOST_CHECK(st == st_other_copy);
}
+
+BOOST_AUTO_TEST_CASE(insert_graph) {
+ std::cout << "********************************************************************" << std::endl;
+ std::cout << "INSERT GRAPH" << std::endl;
+ typedef typename boost::adjacency_list<boost::vecS, boost::vecS,
+ boost::undirectedS,
+ boost::property<vertex_filtration_t, double>,
+ boost::property<edge_filtration_t, double>> Graph;
+ Graph g(3);
+ // filtration value 0 everywhere
+ put(Gudhi::vertex_filtration_t(), g, 0, 0);
+ put(Gudhi::vertex_filtration_t(), g, 1, 0);
+ put(Gudhi::vertex_filtration_t(), g, 2, 0);
+ // vertices don't always occur in sorted order
+ add_edge(0, 1, 0, g);
+ add_edge(2, 1, 0, g);
+ add_edge(2, 0, 0, g);
+
+ typedef Simplex_tree<> typeST;
+ typeST st1;
+ st1.insert_graph(g);
+ BOOST_CHECK(st1.num_simplices() == 6);
+
+ // edges can have multiplicity in the graph unless we replace the first vecS with (hash_)setS
+ add_edge(1, 0, 0, g);
+ add_edge(1, 2, 0, g);
+ add_edge(0, 2, 0, g);
+ add_edge(0, 2, 0, g);
+ typeST st2;
+ st2.insert_graph(g);
+ BOOST_CHECK(st2.num_simplices() == 6);
+}
diff --git a/src/common/include/gudhi/distance_functions.h b/src/common/include/gudhi/distance_functions.h
index c556155e..3a5d1fd5 100644
--- a/src/common/include/gudhi/distance_functions.h
+++ b/src/common/include/gudhi/distance_functions.h
@@ -23,6 +23,10 @@
#ifndef DISTANCE_FUNCTIONS_H_
#define DISTANCE_FUNCTIONS_H_
+#include <gudhi/Debug_utils.h>
+
+#include <boost/range/metafunctions.hpp>
+
#include <cmath> // for std::sqrt
#include <type_traits> // for std::decay
#include <iterator> // for std::begin, std::end
@@ -38,20 +42,29 @@ namespace Gudhi {
* have the same dimension. */
class Euclidean_distance {
public:
+ // boost::range_value is not SFINAE-friendly so we cannot use it in the return type
template< typename Point >
- auto operator()(const Point& p1, const Point& p2) const -> typename std::decay<decltype(*std::begin(p1))>::type {
- auto it1 = p1.begin();
- auto it2 = p2.begin();
- typename Point::value_type dist = 0.;
- for (; it1 != p1.end(); ++it1, ++it2) {
- typename Point::value_type tmp = (*it1) - (*it2);
+ typename std::iterator_traits<typename boost::range_iterator<Point>::type>::value_type
+ operator()(const Point& p1, const Point& p2) const {
+ auto it1 = std::begin(p1);
+ auto it2 = std::begin(p2);
+ typedef typename boost::range_value<Point>::type NT;
+ NT dist = 0;
+ for (; it1 != std::end(p1); ++it1, ++it2) {
+ GUDHI_CHECK(it2 != std::end(p2), "inconsistent point dimensions");
+ NT tmp = *it1 - *it2;
dist += tmp*tmp;
}
- return std::sqrt(dist);
+ GUDHI_CHECK(it2 == std::end(p2), "inconsistent point dimensions");
+ using std::sqrt;
+ return sqrt(dist);
}
template< typename T >
- T operator() (const std::pair< T, T >& f, const std::pair< T, T >& s) {
- return sqrt((f.first-s.first)*(f.first-s.first) + (f.second-s.second)*(f.second-s.second));
+ T operator() (const std::pair< T, T >& f, const std::pair< T, T >& s) const {
+ T dx = f.first - s.first;
+ T dy = f.second - s.second;
+ using std::sqrt;
+ return sqrt(dx*dx + dy*dy);
}
};