diff options
Diffstat (limited to 'src/Simplex_tree')
7 files changed, 246 insertions, 207 deletions
diff --git a/src/Simplex_tree/doc/Intro_simplex_tree.h b/src/Simplex_tree/doc/Intro_simplex_tree.h index 769491d9..6b80d1c9 100644 --- a/src/Simplex_tree/doc/Intro_simplex_tree.h +++ b/src/Simplex_tree/doc/Intro_simplex_tree.h @@ -79,7 +79,6 @@ Number of vertices = 10 Number of simplices = 98 \endcode * 1 incidence relations in a complex. It is consequently faster when accessing the boundary of a simplex, but is less * compact and harder to construct from scratch. * - * \copyright GNU General Public License v3. * @} */ diff --git a/src/Simplex_tree/example/cech_complex_cgal_mini_sphere_3d.cpp b/src/Simplex_tree/example/cech_complex_cgal_mini_sphere_3d.cpp index 217e251f..9bd51106 100644 --- a/src/Simplex_tree/example/cech_complex_cgal_mini_sphere_3d.cpp +++ b/src/Simplex_tree/example/cech_complex_cgal_mini_sphere_3d.cpp @@ -1,5 +1,5 @@ -/* This file is part of the Gudhi Library. The Gudhi library - * (Geometric Understanding in Higher Dimensions) is a generic C++ +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ * library for computational topology. * * Author(s): Clément Maria @@ -33,7 +33,7 @@ #include <string> #include <vector> -#include <limits> // infinity +#include <limits> // infinity #include <utility> // for pair #include <map> @@ -50,15 +50,14 @@ using Vertex_handle = Simplex_tree::Vertex_handle; using Simplex_handle = Simplex_tree::Simplex_handle; using Filtration_value = Simplex_tree::Filtration_value; using Siblings = Simplex_tree::Siblings; -using Graph_t = boost::adjacency_list < boost::vecS, boost::vecS, boost::undirectedS -, boost::property < Gudhi::vertex_filtration_t, Filtration_value > -, boost::property < Gudhi::edge_filtration_t, Filtration_value > ->; -using Edge_t = std::pair< Vertex_handle, Vertex_handle >; +using Graph_t = boost::adjacency_list<boost::vecS, boost::vecS, boost::undirectedS, + boost::property<Gudhi::vertex_filtration_t, Filtration_value>, + boost::property<Gudhi::edge_filtration_t, Filtration_value> >; +using Edge_t = std::pair<Vertex_handle, Vertex_handle>; -using Kernel = CGAL::Epick_d< CGAL::Dimension_tag<3> >; +using Kernel = CGAL::Epick_d<CGAL::Dimension_tag<3> >; using Point = Kernel::Point_d; -using Traits = CGAL::Min_sphere_of_points_d_traits_d<Kernel,Filtration_value,3>; +using Traits = CGAL::Min_sphere_of_points_d_traits_d<Kernel, Filtration_value, 3>; using Min_sphere = CGAL::Min_sphere_of_spheres_d<Traits>; using Points_off_reader = Gudhi::Points_off_reader<Point>; @@ -76,7 +75,7 @@ class Cech_blocker { std::cout << vertex << ", "; #endif // DEBUG_TRACES } - Min_sphere ms(points.begin(),points.end()); + Min_sphere ms(points.begin(), points.end()); Filtration_value radius = ms.radius(); #if DEBUG_TRACES std::cout << "] - radius = " << radius << " - returns " << (radius > threshold_) << std::endl; @@ -85,24 +84,20 @@ class Cech_blocker { return (radius > threshold_); } Cech_blocker(Simplex_tree& simplex_tree, Filtration_value threshold, const std::vector<Point>& point_cloud) - : simplex_tree_(simplex_tree), - threshold_(threshold), - point_cloud_(point_cloud) { } + : simplex_tree_(simplex_tree), threshold_(threshold), point_cloud_(point_cloud) {} + private: Simplex_tree simplex_tree_; Filtration_value threshold_; std::vector<Point> point_cloud_; }; -template< typename InputPointRange> -Graph_t compute_proximity_graph(InputPointRange &points, Filtration_value threshold); +template <typename InputPointRange> +Graph_t compute_proximity_graph(InputPointRange& points, Filtration_value threshold); -void program_options(int argc, char * argv[] - , std::string & off_file_points - , Filtration_value & threshold - , int & dim_max); +void program_options(int argc, char* argv[], std::string& off_file_points, Filtration_value& threshold, int& dim_max); -int main(int argc, char * argv[]) { +int main(int argc, char* argv[]) { std::string off_file_points; Filtration_value threshold; int dim_max; @@ -115,7 +110,7 @@ int main(int argc, char * argv[]) { // Compute the proximity graph of the points Graph_t prox_graph = compute_proximity_graph(off_reader.get_point_cloud(), threshold); - //Min_sphere sph1(off_reader.get_point_cloud()[0], off_reader.get_point_cloud()[1], off_reader.get_point_cloud()[2]); + // Min_sphere sph1(off_reader.get_point_cloud()[0], off_reader.get_point_cloud()[1], off_reader.get_point_cloud()[2]); // Construct the Rips complex in a Simplex Tree Simplex_tree st; // insert the proximity graph in the simplex tree @@ -135,7 +130,8 @@ int main(int argc, char * argv[]) { std::cout << "* The complex contains " << st.num_simplices() << " simplices - dimension=" << st.dimension() << "\n"; std::cout << "* Iterator on Simplices in the filtration, with [filtration value]:\n"; for (auto f_simplex : st.filtration_simplex_range()) { - std::cout << " " << "[" << st.filtration(f_simplex) << "] "; + std::cout << " " + << "[" << st.filtration(f_simplex) << "] "; for (auto vertex : st.simplex_vertex_range(f_simplex)) { std::cout << static_cast<int>(vertex) << " "; } @@ -145,24 +141,19 @@ int main(int argc, char * argv[]) { return 0; } -void program_options(int argc, char * argv[] - , std::string & off_file_points - , Filtration_value & threshold - , int & dim_max) { +void program_options(int argc, char* argv[], std::string& off_file_points, Filtration_value& threshold, int& dim_max) { namespace po = boost::program_options; po::options_description hidden("Hidden options"); - hidden.add_options() - ("input-file", po::value<std::string>(&off_file_points), - "Name of an OFF file containing a 3d point set.\n"); + hidden.add_options()("input-file", po::value<std::string>(&off_file_points), + "Name of an OFF file containing a 3d point set.\n"); po::options_description visible("Allowed options", 100); - visible.add_options() - ("help,h", "produce help message") - ("max-edge-length,r", - po::value<Filtration_value>(&threshold)->default_value(std::numeric_limits<Filtration_value>::infinity()), - "Maximal length of an edge for the Cech complex construction.") - ("cpx-dimension,d", po::value<int>(&dim_max)->default_value(1), - "Maximal dimension of the Cech complex we want to compute."); + visible.add_options()("help,h", "produce help message")( + "max-edge-length,r", + po::value<Filtration_value>(&threshold)->default_value(std::numeric_limits<Filtration_value>::infinity()), + "Maximal length of an edge for the Cech complex construction.")( + "cpx-dimension,d", po::value<int>(&dim_max)->default_value(1), + "Maximal dimension of the Cech complex we want to compute."); po::positional_options_description pos; pos.add("input-file", 1); @@ -171,8 +162,7 @@ void program_options(int argc, char * argv[] all.add(visible).add(hidden); po::variables_map vm; - po::store(po::command_line_parser(argc, argv). - options(all).positional(pos).run(), vm); + po::store(po::command_line_parser(argc, argv).options(all).positional(pos).run(), vm); po::notify(vm); if (vm.count("help") || !vm.count("input-file")) { @@ -194,10 +184,10 @@ void program_options(int argc, char * argv[] * The type PointCloud furnishes .begin() and .end() methods, that return * iterators with value_type Point. */ -template< typename InputPointRange> -Graph_t compute_proximity_graph(InputPointRange &points, Filtration_value threshold) { - std::vector< Edge_t > edges; - std::vector< Filtration_value > edges_fil; +template <typename InputPointRange> +Graph_t compute_proximity_graph(InputPointRange& points, Filtration_value threshold) { + std::vector<Edge_t> edges; + std::vector<Filtration_value> edges_fil; Kernel k; Vertex_handle idx_u, idx_v; @@ -217,16 +207,13 @@ Graph_t compute_proximity_graph(InputPointRange &points, Filtration_value thresh ++idx_u; } - Graph_t skel_graph(edges.begin() - , edges.end() - , edges_fil.begin() - , idx_u); // number of points labeled from 0 to idx_u-1 + Graph_t skel_graph(edges.begin(), edges.end(), edges_fil.begin(), + idx_u); // number of points labeled from 0 to idx_u-1 auto vertex_prop = boost::get(Gudhi::vertex_filtration_t(), skel_graph); boost::graph_traits<Graph_t>::vertex_iterator vi, vi_end; - for (std::tie(vi, vi_end) = boost::vertices(skel_graph); - vi != vi_end; ++vi) { + for (std::tie(vi, vi_end) = boost::vertices(skel_graph); vi != vi_end; ++vi) { boost::put(vertex_prop, *vi, 0.); } diff --git a/src/Simplex_tree/example/graph_expansion_with_blocker.cpp b/src/Simplex_tree/example/graph_expansion_with_blocker.cpp index 86bfb8cb..0d458cbd 100644 --- a/src/Simplex_tree/example/graph_expansion_with_blocker.cpp +++ b/src/Simplex_tree/example/graph_expansion_with_blocker.cpp @@ -27,8 +27,7 @@ using Simplex_tree = Gudhi::Simplex_tree<>; using Simplex_handle = Simplex_tree::Simplex_handle; -int main(int argc, char * const argv[]) { - +int main(int argc, char* const argv[]) { // Construct the Simplex Tree with a 1-skeleton graph example Simplex_tree simplexTree; @@ -45,33 +44,32 @@ int main(int argc, char * const argv[]) { simplexTree.insert_simplex({5, 6}, 10.); simplexTree.insert_simplex({6}, 10.); - simplexTree.expansion_with_blockers(3, [&](Simplex_handle sh){ - bool result = false; - std::cout << "Blocker on ["; - // User can loop on the vertices from the given simplex_handle i.e. - for (auto vertex : simplexTree.simplex_vertex_range(sh)) { - // We block the expansion, if the vertex '6' is in the given list of vertices - if (vertex == 6) - result = true; - std::cout << vertex << ", "; - } - std::cout << "] ( " << simplexTree.filtration(sh); - // User can re-assign a new filtration value directly in the blocker (default is the maximal value of boudaries) - simplexTree.assign_filtration(sh, simplexTree.filtration(sh) + 1.); + simplexTree.expansion_with_blockers(3, [&](Simplex_handle sh) { + bool result = false; + std::cout << "Blocker on ["; + // User can loop on the vertices from the given simplex_handle i.e. + for (auto vertex : simplexTree.simplex_vertex_range(sh)) { + // We block the expansion, if the vertex '6' is in the given list of vertices + if (vertex == 6) result = true; + std::cout << vertex << ", "; + } + std::cout << "] ( " << simplexTree.filtration(sh); + // User can re-assign a new filtration value directly in the blocker (default is the maximal value of boudaries) + simplexTree.assign_filtration(sh, simplexTree.filtration(sh) + 1.); - std::cout << " + 1. ) = " << result << std::endl; + std::cout << " + 1. ) = " << result << std::endl; - return result; - }); + return result; + }); std::cout << "********************************************************************\n"; std::cout << "* The complex contains " << simplexTree.num_simplices() << " simplices"; std::cout << " - dimension " << simplexTree.dimension() << "\n"; std::cout << "* Iterator on Simplices in the filtration, with [filtration value]:\n"; for (auto f_simplex : simplexTree.filtration_simplex_range()) { - std::cout << " " << "[" << simplexTree.filtration(f_simplex) << "] "; - for (auto vertex : simplexTree.simplex_vertex_range(f_simplex)) - std::cout << "(" << vertex << ")"; + std::cout << " " + << "[" << simplexTree.filtration(f_simplex) << "] "; + for (auto vertex : simplexTree.simplex_vertex_range(f_simplex)) std::cout << "(" << vertex << ")"; std::cout << std::endl; } diff --git a/src/Simplex_tree/example/simple_simplex_tree.cpp b/src/Simplex_tree/example/simple_simplex_tree.cpp index b6b65b88..828977c2 100644 --- a/src/Simplex_tree/example/simple_simplex_tree.cpp +++ b/src/Simplex_tree/example/simple_simplex_tree.cpp @@ -30,10 +30,10 @@ using Simplex_tree = Gudhi::Simplex_tree<>; using Vertex_handle = Simplex_tree::Vertex_handle; using Filtration_value = Simplex_tree::Filtration_value; -using typeVectorVertex = std::vector< Vertex_handle >; -using typePairSimplexBool = std::pair< Simplex_tree::Simplex_handle, bool >; +using typeVectorVertex = std::vector<Vertex_handle>; +using typePairSimplexBool = std::pair<Simplex_tree::Simplex_handle, bool>; -int main(int argc, char * const argv[]) { +int main(int argc, char* const argv[]) { const Filtration_value FIRST_FILTRATION_VALUE = 0.1; const Filtration_value SECOND_FILTRATION_VALUE = 0.2; const Filtration_value THIRD_FILTRATION_VALUE = 0.3; @@ -54,7 +54,7 @@ int main(int argc, char * const argv[]) { // ++ FIRST std::cout << " * INSERT 0" << std::endl; - typeVectorVertex firstSimplexVector = { 0 }; + typeVectorVertex firstSimplexVector = {0}; typePairSimplexBool returnValue = simplexTree.insert_simplex(firstSimplexVector, Filtration_value(FIRST_FILTRATION_VALUE)); @@ -66,9 +66,8 @@ int main(int argc, char * const argv[]) { // ++ SECOND std::cout << " * INSERT 1" << std::endl; - typeVectorVertex secondSimplexVector = { 1 }; - returnValue = - simplexTree.insert_simplex(secondSimplexVector, Filtration_value(FIRST_FILTRATION_VALUE)); + typeVectorVertex secondSimplexVector = {1}; + returnValue = simplexTree.insert_simplex(secondSimplexVector, Filtration_value(FIRST_FILTRATION_VALUE)); if (returnValue.second == true) { std::cout << " + 1 INSERTED" << std::endl; @@ -78,9 +77,8 @@ int main(int argc, char * const argv[]) { // ++ THIRD std::cout << " * INSERT (0,1)" << std::endl; - typeVectorVertex thirdSimplexVector = { 0, 1 }; - returnValue = - simplexTree.insert_simplex(thirdSimplexVector, Filtration_value(SECOND_FILTRATION_VALUE)); + typeVectorVertex thirdSimplexVector = {0, 1}; + returnValue = simplexTree.insert_simplex(thirdSimplexVector, Filtration_value(SECOND_FILTRATION_VALUE)); if (returnValue.second == true) { std::cout << " + (0,1) INSERTED" << std::endl; @@ -90,9 +88,8 @@ int main(int argc, char * const argv[]) { // ++ FOURTH std::cout << " * INSERT 2" << std::endl; - typeVectorVertex fourthSimplexVector = { 2 }; - returnValue = - simplexTree.insert_simplex(fourthSimplexVector, Filtration_value(FIRST_FILTRATION_VALUE)); + typeVectorVertex fourthSimplexVector = {2}; + returnValue = simplexTree.insert_simplex(fourthSimplexVector, Filtration_value(FIRST_FILTRATION_VALUE)); if (returnValue.second == true) { std::cout << " + 2 INSERTED" << std::endl; @@ -102,9 +99,8 @@ int main(int argc, char * const argv[]) { // ++ FIFTH std::cout << " * INSERT (2,0)" << std::endl; - typeVectorVertex fifthSimplexVector = { 2, 0 }; - returnValue = - simplexTree.insert_simplex(fifthSimplexVector, Filtration_value(SECOND_FILTRATION_VALUE)); + typeVectorVertex fifthSimplexVector = {2, 0}; + returnValue = simplexTree.insert_simplex(fifthSimplexVector, Filtration_value(SECOND_FILTRATION_VALUE)); if (returnValue.second == true) { std::cout << " + (2,0) INSERTED" << std::endl; @@ -114,9 +110,8 @@ int main(int argc, char * const argv[]) { // ++ SIXTH std::cout << " * INSERT (2,1)" << std::endl; - typeVectorVertex sixthSimplexVector = { 2, 1 }; - returnValue = - simplexTree.insert_simplex(sixthSimplexVector, Filtration_value(SECOND_FILTRATION_VALUE)); + typeVectorVertex sixthSimplexVector = {2, 1}; + returnValue = simplexTree.insert_simplex(sixthSimplexVector, Filtration_value(SECOND_FILTRATION_VALUE)); if (returnValue.second == true) { std::cout << " + (2,1) INSERTED" << std::endl; @@ -126,9 +121,8 @@ int main(int argc, char * const argv[]) { // ++ SEVENTH std::cout << " * INSERT (2,1,0)" << std::endl; - typeVectorVertex seventhSimplexVector = { 2, 1, 0 }; - returnValue = - simplexTree.insert_simplex(seventhSimplexVector, Filtration_value(THIRD_FILTRATION_VALUE)); + typeVectorVertex seventhSimplexVector = {2, 1, 0}; + returnValue = simplexTree.insert_simplex(seventhSimplexVector, Filtration_value(THIRD_FILTRATION_VALUE)); if (returnValue.second == true) { std::cout << " + (2,1,0) INSERTED" << std::endl; @@ -138,9 +132,8 @@ int main(int argc, char * const argv[]) { // ++ EIGHTH std::cout << " * INSERT 3" << std::endl; - typeVectorVertex eighthSimplexVector = { 3 }; - returnValue = - simplexTree.insert_simplex(eighthSimplexVector, Filtration_value(FIRST_FILTRATION_VALUE)); + typeVectorVertex eighthSimplexVector = {3}; + returnValue = simplexTree.insert_simplex(eighthSimplexVector, Filtration_value(FIRST_FILTRATION_VALUE)); if (returnValue.second == true) { std::cout << " + 3 INSERTED" << std::endl; @@ -150,9 +143,8 @@ int main(int argc, char * const argv[]) { // ++ NINETH std::cout << " * INSERT (3,0)" << std::endl; - typeVectorVertex ninethSimplexVector = { 3, 0 }; - returnValue = - simplexTree.insert_simplex(ninethSimplexVector, Filtration_value(SECOND_FILTRATION_VALUE)); + typeVectorVertex ninethSimplexVector = {3, 0}; + returnValue = simplexTree.insert_simplex(ninethSimplexVector, Filtration_value(SECOND_FILTRATION_VALUE)); if (returnValue.second == true) { std::cout << " + (3,0) INSERTED" << std::endl; @@ -162,7 +154,7 @@ int main(int argc, char * const argv[]) { // ++ TENTH std::cout << " * INSERT 0 (already inserted)" << std::endl; - typeVectorVertex tenthSimplexVector = { 0 }; + typeVectorVertex tenthSimplexVector = {0}; // With a different filtration value returnValue = simplexTree.insert_simplex(tenthSimplexVector, Filtration_value(FOURTH_FILTRATION_VALUE)); @@ -174,9 +166,8 @@ int main(int argc, char * const argv[]) { // ++ ELEVENTH std::cout << " * INSERT (2,1,0) (already inserted)" << std::endl; - typeVectorVertex eleventhSimplexVector = { 2, 1, 0 }; - returnValue = - simplexTree.insert_simplex(eleventhSimplexVector, Filtration_value(FOURTH_FILTRATION_VALUE)); + typeVectorVertex eleventhSimplexVector = {2, 1, 0}; + returnValue = simplexTree.insert_simplex(eleventhSimplexVector, Filtration_value(FOURTH_FILTRATION_VALUE)); if (returnValue.second == true) { std::cout << " + (2,1,0) INSERTED" << std::endl; @@ -192,9 +183,9 @@ int main(int argc, char * const argv[]) { std::cout << " - dimension " << simplexTree.dimension() << "\n"; std::cout << "* Iterator on Simplices in the filtration, with [filtration value]:\n"; for (auto f_simplex : simplexTree.filtration_simplex_range()) { - std::cout << " " << "[" << simplexTree.filtration(f_simplex) << "] "; - for (auto vertex : simplexTree.simplex_vertex_range(f_simplex)) - std::cout << "(" << vertex << ")"; + std::cout << " " + << "[" << simplexTree.filtration(f_simplex) << "] "; + for (auto vertex : simplexTree.simplex_vertex_range(f_simplex)) std::cout << "(" << vertex << ")"; std::cout << std::endl; } // [0.1] 0 @@ -217,7 +208,7 @@ int main(int argc, char * const argv[]) { else std::cout << "***- NO IT ISN'T\n"; - typeVectorVertex unknownSimplexVector = { 15 }; + typeVectorVertex unknownSimplexVector = {15}; simplexFound = simplexTree.find(unknownSimplexVector); std::cout << "**************IS THE SIMPLEX {15} IN THE SIMPLEX TREE ?\n"; if (simplexFound != simplexTree.null_simplex()) @@ -232,7 +223,7 @@ int main(int argc, char * const argv[]) { else std::cout << "***- NO IT ISN'T\n"; - typeVectorVertex otherSimplexVector = { 1, 15 }; + typeVectorVertex otherSimplexVector = {1, 15}; simplexFound = simplexTree.find(otherSimplexVector); std::cout << "**************IS THE SIMPLEX {15,1} IN THE SIMPLEX TREE ?\n"; if (simplexFound != simplexTree.null_simplex()) @@ -240,7 +231,7 @@ int main(int argc, char * const argv[]) { else std::cout << "***- NO IT ISN'T\n"; - typeVectorVertex invSimplexVector = { 1, 2, 0 }; + typeVectorVertex invSimplexVector = {1, 2, 0}; simplexFound = simplexTree.find(invSimplexVector); std::cout << "**************IS THE SIMPLEX {1,2,0} IN THE SIMPLEX TREE ?\n"; if (simplexFound != simplexTree.null_simplex()) @@ -248,7 +239,7 @@ int main(int argc, char * const argv[]) { else std::cout << "***- NO IT ISN'T\n"; - simplexFound = simplexTree.find({ 0, 1 }); + simplexFound = simplexTree.find({0, 1}); std::cout << "**************IS THE SIMPLEX {0,1} IN THE SIMPLEX TREE ?\n"; if (simplexFound != simplexTree.null_simplex()) std::cout << "***+ YES IT IS!\n"; @@ -256,23 +247,20 @@ int main(int argc, char * const argv[]) { std::cout << "***- NO IT ISN'T\n"; std::cout << "**************COFACES OF {0,1} IN CODIMENSION 1 ARE\n"; - for (auto& simplex : simplexTree.cofaces_simplex_range(simplexTree.find({0,1}), 1)) { - for (auto vertex : simplexTree.simplex_vertex_range(simplex)) - std::cout << "(" << vertex << ")"; + for (auto& simplex : simplexTree.cofaces_simplex_range(simplexTree.find({0, 1}), 1)) { + for (auto vertex : simplexTree.simplex_vertex_range(simplex)) std::cout << "(" << vertex << ")"; std::cout << std::endl; } std::cout << "**************STARS OF {0,1} ARE\n"; - for (auto& simplex : simplexTree.star_simplex_range(simplexTree.find({0,1}))) { - for (auto vertex : simplexTree.simplex_vertex_range(simplex)) - std::cout << "(" << vertex << ")"; + for (auto& simplex : simplexTree.star_simplex_range(simplexTree.find({0, 1}))) { + for (auto vertex : simplexTree.simplex_vertex_range(simplex)) std::cout << "(" << vertex << ")"; std::cout << std::endl; } std::cout << "**************BOUNDARIES OF {0,1,2} ARE\n"; - for (auto& simplex : simplexTree.boundary_simplex_range(simplexTree.find({0,1,2}))) { - for (auto vertex : simplexTree.simplex_vertex_range(simplex)) - std::cout << "(" << vertex << ")"; + for (auto& simplex : simplexTree.boundary_simplex_range(simplexTree.find({0, 1, 2}))) { + for (auto vertex : simplexTree.simplex_vertex_range(simplex)) std::cout << "(" << vertex << ")"; std::cout << std::endl; } 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; diff --git a/src/Simplex_tree/include/gudhi/Simplex_tree/Simplex_tree_iterators.h b/src/Simplex_tree/include/gudhi/Simplex_tree/Simplex_tree_iterators.h index 7e0a454d..ab7346d4 100644 --- a/src/Simplex_tree/include/gudhi/Simplex_tree/Simplex_tree_iterators.h +++ b/src/Simplex_tree/include/gudhi/Simplex_tree/Simplex_tree_iterators.h @@ -23,6 +23,8 @@ #ifndef SIMPLEX_TREE_SIMPLEX_TREE_ITERATORS_H_ #define SIMPLEX_TREE_SIMPLEX_TREE_ITERATORS_H_ +#include <gudhi/Debug_utils.h> + #include <boost/iterator/iterator_facade.hpp> #include <boost/version.hpp> #if BOOST_VERSION >= 105600 @@ -109,11 +111,18 @@ class Simplex_tree_boundary_simplex_iterator : public boost::iterator_facade< : last_(sh->first), sib_(nullptr), st_(st) { + // Only check once at the beginning instead of for every increment, as this is expensive. + if (SimplexTree::Options::contiguous_vertices) + GUDHI_CHECK(st_->contiguous_vertices(), "The set of vertices is not { 0, ..., n } without holes"); Siblings * sib = st->self_siblings(sh); next_ = sib->parent(); sib_ = sib->oncles(); if (sib_ != nullptr) { - sh_ = sib_->find(next_); + if (SimplexTree::Options::contiguous_vertices && sib_->oncles() == nullptr) + // Only relevant for edges + sh_ = sib_->members_.begin()+next_; + else + sh_ = sib_->find(next_); } else { sh_ = st->null_simplex(); } // vertex: == end() @@ -140,14 +149,19 @@ class Simplex_tree_boundary_simplex_iterator : public boost::iterator_facade< Siblings * for_sib = sib_; Siblings * new_sib = sib_->oncles(); auto rit = suffix_.rbegin(); - if (SimplexTree::Options::contiguous_vertices && new_sib == nullptr && rit != suffix_.rend()) { - // We reached the root, use a short-cut to find a vertex. We could also - // optimize finding the second vertex of a segment, but people are - // expected to call endpoints(). - assert(st_->contiguous_vertices()); - sh_ = for_sib->members_.begin()+*rit; - for_sib = sh_->second.children(); - ++rit; + if (SimplexTree::Options::contiguous_vertices && new_sib == nullptr) { + // We reached the root, use a short-cut to find a vertex. + if (rit == suffix_.rend()) { + // Segment, this vertex is the last boundary simplex + sh_ = for_sib->members_.begin()+last_; + sib_ = nullptr; + return; + } else { + // Dim >= 2, initial step of the descent + sh_ = for_sib->members_.begin()+*rit; + for_sib = sh_->second.children(); + ++rit; + } } for (; rit != suffix_.rend(); ++rit) { sh_ = for_sib->find(*rit); 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); +} |