From fa1ff1284f51341d452b6adc3a667d13a5d35747 Mon Sep 17 00:00:00 2001 From: fgodi Date: Mon, 4 Apr 2016 13:14:31 +0000 Subject: cmake lists git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/bottleneckDistance@1092 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: 77b2a427cb3ff80b9a0c5d1c71f9f1de11547530 --- CMakeLists.txt | 19 +++++++++++++++---- 1 file changed, 15 insertions(+), 4 deletions(-) (limited to 'CMakeLists.txt') diff --git a/CMakeLists.txt b/CMakeLists.txt index f30c1785..6ed0aedd 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -42,7 +42,6 @@ if (GPROF_PATH) message("gprof found in ${GPROF_PATH}") endif() - if(NOT Boost_FOUND) message(FATAL_ERROR "NOTICE: This demo requires Boost and will not be compiled.") else() @@ -75,12 +74,24 @@ else() add_subdirectory(src/Hasse_complex/example) add_subdirectory(src/Alpha_shapes/example) add_subdirectory(src/Alpha_shapes/test) - add_subdirectory(src/Bottleneck/example) - add_subdirectory(src/Bottleneck/test) + add_subdirectory(src/Bipartite_graphs_matching/example) + add_subdirectory(src/Bipartite_graphs_matching/test) # data points generator add_subdirectory(data/points/generator) -endif() +endif() +include( ${CGAL_USE_FILE} ) +# In CMakeLists.txt, when include(${CGAL_USE_FILE}), CXX_FLAGS are overwritten. +# cf. http://doc.cgal.org/latest/Manual/installation.html#title40 +# A workaround is to add "-std=c++11" again. +if(NOT MSVC) +include(CheckCXXCompilerFlag) +CHECK_CXX_COMPILER_FLAG(-std=c++11 COMPILER_SUPPORTS_CXX11) +if(COMPILER_SUPPORTS_CXX11) +set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++11") +endif() +endif() +# - End of workaround -- cgit v1.2.3 From 10ee730fd273592752d3e05cf8805f9d68efa5b4 Mon Sep 17 00:00:00 2001 From: fgodi Date: Mon, 4 Apr 2016 13:57:11 +0000 Subject: Strange : compile only without the workaround c++11, then removed git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/bottleneckDistance@1095 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: d1fd8df7a0aeccd69f029b8848c117cb3b567bd1 --- CMakeLists.txt | 13 ------------- 1 file changed, 13 deletions(-) (limited to 'CMakeLists.txt') diff --git a/CMakeLists.txt b/CMakeLists.txt index 6ed0aedd..54165c59 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -82,16 +82,3 @@ else() endif() -include( ${CGAL_USE_FILE} ) -# In CMakeLists.txt, when include(${CGAL_USE_FILE}), CXX_FLAGS are overwritten. -# cf. http://doc.cgal.org/latest/Manual/installation.html#title40 -# A workaround is to add "-std=c++11" again. -if(NOT MSVC) -include(CheckCXXCompilerFlag) -CHECK_CXX_COMPILER_FLAG(-std=c++11 COMPILER_SUPPORTS_CXX11) -if(COMPILER_SUPPORTS_CXX11) -set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++11") -endif() -endif() -# - End of workaround - -- cgit v1.2.3 From 9db309ff6abbcd2ba3ebbc74fa387b55d9660789 Mon Sep 17 00:00:00 2001 From: fgodi Date: Wed, 25 May 2016 14:02:14 +0000 Subject: save git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/bottleneckDistance@1199 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: 00868265725afaf47ed8292c0d6f7cff235a016f --- CMakeLists.txt | 2 +- src/Bipartite_graphs_matching/example/basic.cpp | 42 +++++++++++++++++++++- .../include/gudhi/Planar_neighbors_finder.h | 11 +++--- .../test/bottleneck_unit_test.cpp | 37 ------------------- 4 files changed, 46 insertions(+), 46 deletions(-) (limited to 'CMakeLists.txt') diff --git a/CMakeLists.txt b/CMakeLists.txt index 54165c59..0f031a49 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -14,7 +14,7 @@ if(MSVC) # Turn off some VC++ warnings SET (CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} /wd4267 /wd4668 /wd4311 /wd4800 /wd4820 /wd4503 /wd4244 /wd4345 /wd4996 /wd4396 /wd4018") else() - set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -O2 -std=c++11 -Wall -Wpedantic -Wsign-compare") + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -O0 -g -std=c++11 -Wall -Wpedantic -Wsign-compare") set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -ggdb -O0") set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE}") endif() diff --git a/src/Bipartite_graphs_matching/example/basic.cpp b/src/Bipartite_graphs_matching/example/basic.cpp index 772bcd4a..797968cc 100644 --- a/src/Bipartite_graphs_matching/example/basic.cpp +++ b/src/Bipartite_graphs_matching/example/basic.cpp @@ -23,8 +23,48 @@ #include "../include/gudhi/Graph_matching.h" #include +#include +#include + using namespace Gudhi::bipartite_graph_matching; + +double upper_bound = 400.; // any real >0 + +int main(){ + std::ofstream objetfichier; + objetfichier.open("results.csv", std::ios::out); + + for(int n =50; n<=1000; n+=100){ + std::uniform_real_distribution unif1(0.,upper_bound); + std::uniform_real_distribution unif2(upper_bound/1000.,upper_bound/100.); + std::default_random_engine re; + std::vector< std::pair > v1, v2; + for (int i = 0; i < n; i++) { + double a = unif1(re); + double b = unif1(re); + double x = unif2(re); + double y = unif2(re); + v1.emplace_back(std::min(a,b), std::max(a,b)); + v2.emplace_back(std::min(a,b)+std::min(x,y), std::max(a,b)+std::max(x,y)); + if(i%5==0) + v1.emplace_back(std::min(a,b),std::min(a,b)+x); + if(i%3==0) + v2.emplace_back(std::max(a,b),std::max(a,b)+y); + } + + std::chrono::steady_clock::time_point start = std::chrono::steady_clock::now(); + double b = bottleneck_distance(v1,v2); + std::chrono::steady_clock::time_point end = std::chrono::steady_clock::now(); + + typedef std::chrono::duration millisecs_t; + millisecs_t duration(std::chrono::duration_cast(end-start)); + objetfichier << n << ";" << duration.count() << ";" << b << std::endl; + } + objetfichier.close(); +} + +/* int main() { std::vector< std::pair > v1, v2; @@ -40,4 +80,4 @@ int main() { std::cout << "Bottleneck distance = " << b << std::endl; -} +}*/ diff --git a/src/Bipartite_graphs_matching/include/gudhi/Planar_neighbors_finder.h b/src/Bipartite_graphs_matching/include/gudhi/Planar_neighbors_finder.h index abb14e29..8231409c 100644 --- a/src/Bipartite_graphs_matching/include/gudhi/Planar_neighbors_finder.h +++ b/src/Bipartite_graphs_matching/include/gudhi/Planar_neighbors_finder.h @@ -82,7 +82,7 @@ public: void remove(int v_point_index); /** \internal \brief Can the point given as parameter be returned ? */ bool contains(int v_point_index) const; - /** \internal \brief Provide and remove a V point near to the U point given as parameter, null_point_index() if there isn't such a point. */ + /** \internal \brief Provide a V point near to the U point given as parameter, null_point_index() if there isn't such a point. */ int pull_near(int u_point_index); /** \internal \brief Provide and remove all the V points near to the U point given as parameter. */ virtual std::shared_ptr< std::list > pull_all_near(int u_point_index); @@ -110,6 +110,7 @@ inline void Naive_pnf::add(int v_point_index) { } inline void Naive_pnf::remove(int v_point_index) { + if(!contains(v_point_index)) for(auto it = grid.find(get_v_key(v_point_index)); it!=grid.end(); it++) if(it->second==v_point_index){ grid.erase(it); @@ -135,7 +136,6 @@ inline int Naive_pnf::pull_near(int u_point_index) { for(auto it = grid.find(std::make_pair(i0 +(i%3)-1, j0+(j%3)-1)); it!=grid.end(); it++) if (G::distance(u_point_index, it->second) <= r) { int tmp = it->second; - grid.erase(it); return tmp; } return null_point_index(); @@ -186,12 +186,8 @@ inline int Cgal_pnf::pull_near(int u_point_index){ std::vector w = {1., 1.}; K_neighbor_search search(kd_t, u_point, 0., true, Distance(0, 2, w)); auto it = search.begin(); - if(it==search.end() || G::distance(u_point_index, it->first.point_index) > r){ - for(auto itc=contents.cbegin(); itc != contents.cend(); itc++) - if(G::distance(u_point_index, *itc) <= r) - std::cout << G::distance(u_point_index, *itc) << " ! > " << r << std::endl; + if(it==search.end() || G::distance(u_point_index, it->first.point_index) > r) return null_point_index(); - } return it->first.point_index; } @@ -200,6 +196,7 @@ inline std::shared_ptr< std::list > Cgal_pnf::pull_all_near(int u_point_ind int last_pull = pull_near(u_point_index); while (last_pull != null_point_index()) { all_pull->emplace_back(last_pull); + remove(last_pull); last_pull = pull_near(u_point_index); } return all_pull; diff --git a/src/Bipartite_graphs_matching/test/bottleneck_unit_test.cpp b/src/Bipartite_graphs_matching/test/bottleneck_unit_test.cpp index d53a5a3c..7c1d5614 100644 --- a/src/Bipartite_graphs_matching/test/bottleneck_unit_test.cpp +++ b/src/Bipartite_graphs_matching/test/bottleneck_unit_test.cpp @@ -4,9 +4,6 @@ #include #include "../include/gudhi/Graph_matching.h" -#include -#include - using namespace Gudhi::bipartite_graph_matching; int n1 = 81; // a natural number >0 @@ -166,37 +163,3 @@ BOOST_AUTO_TEST_CASE(global){ } BOOST_CHECK(bottleneck_distance(v1, v2) <= upper_bound/100.); } - -/* -BOOST_AUTO_TEST_CASE(chrono) { - std::ofstream objetfichier; - objetfichier.open("results.csv", std::ios::out); - - for(int n =50; n<=1000; n+=100){ - std::uniform_real_distribution unif1(0.,upper_bound); - std::uniform_real_distribution unif2(upper_bound/1000.,upper_bound/100.); - std::default_random_engine re; - std::vector< std::pair > v1, v2; - for (int i = 0; i < n; i++) { - double a = unif1(re); - double b = unif1(re); - double x = unif2(re); - double y = unif2(re); - v1.emplace_back(std::min(a,b), std::max(a,b)); - v2.emplace_back(std::min(a,b)+std::min(x,y), std::max(a,b)+std::max(x,y)); - if(i%5==0) - v1.emplace_back(std::min(a,b),std::min(a,b)+x); - if(i%3==0) - v2.emplace_back(std::max(a,b),std::max(a,b)+y); - } - - std::chrono::steady_clock::time_point start = std::chrono::steady_clock::now(); - double b = bottleneck_distance(v1,v2); - std::chrono::steady_clock::time_point end = std::chrono::steady_clock::now(); - - typedef std::chrono::duration millisecs_t; - millisecs_t duration(std::chrono::duration_cast(end-start)); - objetfichier << n << ";" << duration.count() << ";" << b << std::endl; - } - objetfichier.close(); -}*/ -- cgit v1.2.3 From 40a293e774742d308ce02416326e67d5375cc371 Mon Sep 17 00:00:00 2001 From: cjamin Date: Wed, 25 May 2016 15:31:07 +0000 Subject: Missing include_directories git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/bottleneckDistance@1203 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: f75268759affc6ce3e7c5a2929ba00971a1c48a1 --- CMakeLists.txt | 1 + src/Bipartite_graphs_matching/example/basic.cpp | 2 +- 2 files changed, 2 insertions(+), 1 deletion(-) (limited to 'CMakeLists.txt') diff --git a/CMakeLists.txt b/CMakeLists.txt index 0f031a49..066a49cf 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -57,6 +57,7 @@ else() include_directories(src/common/include/) include_directories(src/Alpha_shapes/include/) + include_directories(src/Bipartite_graphs_matching/include/) include_directories(src/Bottleneck/include/) include_directories(src/Contraction/include/) include_directories(src/Hasse_complex/include/) diff --git a/src/Bipartite_graphs_matching/example/basic.cpp b/src/Bipartite_graphs_matching/example/basic.cpp index 797968cc..f8585fe6 100644 --- a/src/Bipartite_graphs_matching/example/basic.cpp +++ b/src/Bipartite_graphs_matching/example/basic.cpp @@ -20,7 +20,7 @@ * along with this program. If not, see . */ -#include "../include/gudhi/Graph_matching.h" +#include #include #include -- cgit v1.2.3 From bd9b3b687b715a968216032f20f08fbaea7f3baa Mon Sep 17 00:00:00 2001 From: fgodi Date: Thu, 2 Jun 2016 17:16:04 +0000 Subject: lists git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/bottleneckDistance@1242 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: 134d7d5c1f940ac342eb1721de6ee9cb41f1102b --- CMakeLists.txt | 6 +++--- src/Bottleneck_distance/example/CMakeLists.txt | 2 +- src/Bottleneck_distance/test/CMakeLists.txt | 3 ++- 3 files changed, 6 insertions(+), 5 deletions(-) (limited to 'CMakeLists.txt') diff --git a/CMakeLists.txt b/CMakeLists.txt index 066a49cf..dccd59d8 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -57,7 +57,7 @@ else() include_directories(src/common/include/) include_directories(src/Alpha_shapes/include/) - include_directories(src/Bipartite_graphs_matching/include/) + include_directories(src/Bottleneck_distance/include/) include_directories(src/Bottleneck/include/) include_directories(src/Contraction/include/) include_directories(src/Hasse_complex/include/) @@ -75,8 +75,8 @@ else() add_subdirectory(src/Hasse_complex/example) add_subdirectory(src/Alpha_shapes/example) add_subdirectory(src/Alpha_shapes/test) - add_subdirectory(src/Bipartite_graphs_matching/example) - add_subdirectory(src/Bipartite_graphs_matching/test) + add_subdirectory(src/Bottleneck_distance/example) + add_subdirectory(src/Bottleneck_distance/test) # data points generator add_subdirectory(data/points/generator) diff --git a/src/Bottleneck_distance/example/CMakeLists.txt b/src/Bottleneck_distance/example/CMakeLists.txt index 4e8766cc..cc1b9091 100644 --- a/src/Bottleneck_distance/example/CMakeLists.txt +++ b/src/Bottleneck_distance/example/CMakeLists.txt @@ -24,7 +24,7 @@ if(CGAL_FOUND) include( ${EIGEN3_USE_FILE} ) include_directories (BEFORE "../../include") - add_executable (basic_example basic.cpp) + add_executable (bottleneck_example bottleneck_example.cpp) #add_test(dtoffrw_tore3D ${CMAKE_CURRENT_BINARY_DIR}/dtoffrw ${CMAKE_SOURCE_DIR}/data/points/tore3D_1307.off 3) else() diff --git a/src/Bottleneck_distance/test/CMakeLists.txt b/src/Bottleneck_distance/test/CMakeLists.txt index ab0461cc..1338ed6b 100644 --- a/src/Bottleneck_distance/test/CMakeLists.txt +++ b/src/Bottleneck_distance/test/CMakeLists.txt @@ -44,7 +44,8 @@ endif() set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -ggdb -O0") set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE}") -add_executable ( BottleneckUT bottleneck_unit_test.cpp ) +add_executable ( bottleneckUT bottleneck_unit_test.cpp ) +add_executable ( bottleneck_chrono bottleneck_chrono.cpp ) target_link_libraries(BottleneckUT ${Boost_SYSTEM_LIBRARY} ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY}) # Unitary tests -- cgit v1.2.3 From ca6b3d72f2114cdd2d0899fd44dba19303dd3bb2 Mon Sep 17 00:00:00 2001 From: vrouvrea Date: Wed, 28 Sep 2016 11:15:34 +0000 Subject: Add rips complex module and first modification step git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/rips_complex_module@1577 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: fc949570c6cffa4756f84a15c41fadc8c45b2af2 --- CMakeLists.txt | 2 + .../example/rips_multifield_persistence.cpp | 78 +++-- .../example/rips_persistence.cpp | 79 +++-- .../concept/Simplicial_complex_for_rips.h | 53 ++++ src/Rips_complex/doc/Intro_rips_complex.h | 104 +++++++ .../doc/rips_complex_representation.ipe | 326 +++++++++++++++++++++ .../doc/rips_complex_representation.png | Bin 0 -> 15677 bytes src/Rips_complex/doc/rips_one_skeleton.ipe | 326 +++++++++++++++++++++ src/Rips_complex/doc/rips_one_skeleton.png | Bin 0 -> 47651 bytes src/Rips_complex/example/CMakeLists.txt | 8 + .../example/example_rips_complex_from_off_file.cpp | 71 +++++ src/Rips_complex/include/gudhi/Rips_complex.h | 154 ++++++++++ src/Simplex_tree/include/gudhi/Simplex_tree.h | 2 +- src/common/include/gudhi/distance_functions.h | 2 +- .../include/gudhi/graph_simplicial_complex.h | 8 +- 15 files changed, 1117 insertions(+), 96 deletions(-) create mode 100644 src/Rips_complex/concept/Simplicial_complex_for_rips.h create mode 100644 src/Rips_complex/doc/Intro_rips_complex.h create mode 100644 src/Rips_complex/doc/rips_complex_representation.ipe create mode 100644 src/Rips_complex/doc/rips_complex_representation.png create mode 100644 src/Rips_complex/doc/rips_one_skeleton.ipe create mode 100644 src/Rips_complex/doc/rips_one_skeleton.png create mode 100644 src/Rips_complex/example/CMakeLists.txt create mode 100644 src/Rips_complex/example/example_rips_complex_from_off_file.cpp create mode 100644 src/Rips_complex/include/gudhi/Rips_complex.h (limited to 'CMakeLists.txt') diff --git a/CMakeLists.txt b/CMakeLists.txt index 42a56fdc..89e130f4 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -123,6 +123,7 @@ else() include_directories(src/Contraction/include/) include_directories(src/Hasse_complex/include/) include_directories(src/Persistent_cohomology/include/) + include_directories(src/Rips_complex/include/) include_directories(src/Simplex_tree/include/) include_directories(src/Skeleton_blocker/include/) include_directories(src/Witness_complex/include/) @@ -142,6 +143,7 @@ else() add_subdirectory(src/Bitmap_cubical_complex/example) add_subdirectory(src/Alpha_complex/example) add_subdirectory(src/Alpha_complex/test) + add_subdirectory(src/Rips_complex/example) # data points generator add_subdirectory(data/points/generator) diff --git a/src/Persistent_cohomology/example/rips_multifield_persistence.cpp b/src/Persistent_cohomology/example/rips_multifield_persistence.cpp index c5cd775d..d4a07bb6 100644 --- a/src/Persistent_cohomology/example/rips_multifield_persistence.cpp +++ b/src/Persistent_cohomology/example/rips_multifield_persistence.cpp @@ -20,8 +20,7 @@ * along with this program. If not, see . */ -#include -#include +#include #include #include #include @@ -32,14 +31,11 @@ #include #include -using namespace Gudhi; -using namespace Gudhi::persistent_cohomology; -typedef int Vertex_handle; typedef double Filtration_value; void program_options(int argc, char * argv[] - , std::string & filepoints + , std::string & off_file_points , std::string & filediag , Filtration_value & threshold , int & dim_max @@ -48,7 +44,7 @@ void program_options(int argc, char * argv[] , Filtration_value & min_persistence); int main(int argc, char * argv[]) { - std::string filepoints; + std::string off_file_points; std::string filediag; Filtration_value threshold; int dim_max; @@ -56,49 +52,43 @@ int main(int argc, char * argv[]) { int max_p; Filtration_value min_persistence; - program_options(argc, argv, filepoints, filediag, threshold, dim_max, min_p, max_p, min_persistence); + program_options(argc, argv, off_file_points, filediag, threshold, dim_max, min_p, max_p, min_persistence); - // Extract the points from the file filepoints - typedef std::vector Point_t; - std::vector< Point_t > points; - read_points(filepoints, points); - - // Compute the proximity graph of the points - Graph_t prox_graph = compute_proximity_graph(points, threshold - , euclidean_distance); + Gudhi::rips_complex::Rips_complex<> rips_complex_from_file(off_file_points); // Construct the Rips complex in a Simplex Tree - typedef Simplex_tree ST; - ST st; - // insert the proximity graph in the simplex tree - st.insert_graph(prox_graph); - // expand the graph until dimension dim_max - st.expansion(dim_max); - - // Sort the simplices in the order of the filtration - st.initialize_filtration(); - - // Compute the persistence diagram of the complex - Persistent_cohomology pcoh(st); - // initializes the coefficient field for homology - pcoh.init_coefficients(min_p, max_p); - // compute persistent homology, disgarding persistent features of life shorter than min_persistence - pcoh.compute_persistent_cohomology(min_persistence); - - // Output the diagram in filediag - if (filediag.empty()) { - pcoh.output_diagram(); - } else { - std::ofstream out(filediag); - pcoh.output_diagram(out); - out.close(); + using Simplex_tree = Gudhi::Simplex_tree; + Simplex_tree simplex_tree; + + if (rips_complex_from_file.create_complex(simplex_tree, threshold, dim_max, euclidean_distance>)) { + std::cout << "The complex contains " << simplex_tree.num_simplices() << " simplices \n"; + std::cout << " and has dimension " << simplex_tree.dimension() << " \n"; + + // Sort the simplices in the order of the filtration + simplex_tree.initialize_filtration(); + + using Multi_field = Gudhi::persistent_cohomology::Multi_field; + // Compute the persistence diagram of the complex + Gudhi::persistent_cohomology::Persistent_cohomology pcoh(simplex_tree); + // initializes the coefficient field for homology + pcoh.init_coefficients(min_p, max_p); + + pcoh.compute_persistent_cohomology(min_persistence); + + // Output the diagram in filediag + if (filediag.empty()) { + pcoh.output_diagram(); + } else { + std::ofstream out(filediag); + pcoh.output_diagram(out); + out.close(); + } } - return 0; } void program_options(int argc, char * argv[] - , std::string & filepoints + , std::string & off_file_points , std::string & filediag , Filtration_value & threshold , int & dim_max @@ -108,8 +98,8 @@ void program_options(int argc, char * argv[] namespace po = boost::program_options; po::options_description hidden("Hidden options"); hidden.add_options() - ("input-file", po::value(&filepoints), - "Name of file containing a point set. Format is one point per line: X1 ... Xd \n"); + ("input-file", po::value(&off_file_points), + "Name of an OFF file containing a point set.\n"); po::options_description visible("Allowed options"); visible.add_options() diff --git a/src/Persistent_cohomology/example/rips_persistence.cpp b/src/Persistent_cohomology/example/rips_persistence.cpp index cab49395..e3ab5927 100644 --- a/src/Persistent_cohomology/example/rips_persistence.cpp +++ b/src/Persistent_cohomology/example/rips_persistence.cpp @@ -20,7 +20,7 @@ * along with this program. If not, see . */ -#include +#include #include #include #include @@ -32,14 +32,10 @@ #include #include // infinity -using namespace Gudhi; -using namespace Gudhi::persistent_cohomology; - -typedef int Vertex_handle; typedef double Filtration_value; void program_options(int argc, char * argv[] - , std::string & filepoints + , std::string & off_file_points , std::string & filediag , Filtration_value & threshold , int & dim_max @@ -47,59 +43,50 @@ void program_options(int argc, char * argv[] , Filtration_value & min_persistence); int main(int argc, char * argv[]) { - std::string filepoints; + std::string off_file_points; std::string filediag; Filtration_value threshold; int dim_max; int p; Filtration_value min_persistence; - program_options(argc, argv, filepoints, filediag, threshold, dim_max, p, min_persistence); - - // Extract the points from the file filepoints - typedef std::vector Point_t; - std::vector< Point_t > points; - read_points(filepoints, points); + program_options(argc, argv, off_file_points, filediag, threshold, dim_max, p, min_persistence); - // Compute the proximity graph of the points - Graph_t prox_graph = compute_proximity_graph(points, threshold - , euclidean_distance); + Gudhi::rips_complex::Rips_complex<> rips_complex_from_file(off_file_points); // Construct the Rips complex in a Simplex Tree - typedef Simplex_tree ST; - ST st; - // insert the proximity graph in the simplex tree - st.insert_graph(prox_graph); - // expand the graph until dimension dim_max - st.expansion(dim_max); - - std::cout << "The complex contains " << st.num_simplices() << " simplices \n"; - std::cout << " and has dimension " << st.dimension() << " \n"; - - // Sort the simplices in the order of the filtration - st.initialize_filtration(); - - // Compute the persistence diagram of the complex - persistent_cohomology::Persistent_cohomology pcoh(st); - // initializes the coefficient field for homology - pcoh.init_coefficients(p); - - pcoh.compute_persistent_cohomology(min_persistence); - - // Output the diagram in filediag - if (filediag.empty()) { - pcoh.output_diagram(); - } else { - std::ofstream out(filediag); - pcoh.output_diagram(out); - out.close(); + using Simplex_tree = Gudhi::Simplex_tree; + Simplex_tree simplex_tree; + + if (rips_complex_from_file.create_complex(simplex_tree, threshold, dim_max, euclidean_distance>)) { + std::cout << "The complex contains " << simplex_tree.num_simplices() << " simplices \n"; + std::cout << " and has dimension " << simplex_tree.dimension() << " \n"; + + // Sort the simplices in the order of the filtration + simplex_tree.initialize_filtration(); + + using Field_Zp = Gudhi::persistent_cohomology::Field_Zp; + // Compute the persistence diagram of the complex + Gudhi::persistent_cohomology::Persistent_cohomology pcoh(simplex_tree); + // initializes the coefficient field for homology + pcoh.init_coefficients(p); + + pcoh.compute_persistent_cohomology(min_persistence); + + // Output the diagram in filediag + if (filediag.empty()) { + pcoh.output_diagram(); + } else { + std::ofstream out(filediag); + pcoh.output_diagram(out); + out.close(); + } } - return 0; } void program_options(int argc, char * argv[] - , std::string & filepoints + , std::string & off_file_points , std::string & filediag , Filtration_value & threshold , int & dim_max @@ -108,7 +95,7 @@ void program_options(int argc, char * argv[] namespace po = boost::program_options; po::options_description hidden("Hidden options"); hidden.add_options() - ("input-file", po::value(&filepoints), + ("input-file", po::value(&off_file_points), "Name of file containing a point set. Format is one point per line: X1 ... Xd "); po::options_description visible("Allowed options", 100); diff --git a/src/Rips_complex/concept/Simplicial_complex_for_rips.h b/src/Rips_complex/concept/Simplicial_complex_for_rips.h new file mode 100644 index 00000000..470860e9 --- /dev/null +++ b/src/Rips_complex/concept/Simplicial_complex_for_rips.h @@ -0,0 +1,53 @@ +/* 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): Vincent Rouvreau + * + * Copyright (C) 2016 INRIA + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ + +#ifndef CONCEPT_RIPS_COMPLEX_SIMPLICIAL_COMPLEX_FOR_RIPS_H_ +#define CONCEPT_RIPS_COMPLEX_SIMPLICIAL_COMPLEX_FOR_RIPS_H_ + +namespace Gudhi { + +namespace rips_complex { + +/** \brief The concept SimplicialComplexForRips describes the requirements for a type to implement a simplicial + * complex, that can be created from a `Rips_complex`. + */ +struct SimplicialComplexForRips { + /** Handle to specify the simplex filtration value. */ + typedef unspecified Filtration_value; + + /** Returns the number of vertices in the simplicial complex. */ + std::size_t num_vertices(); + + /** \brief Inserts a a given range 'OneSkeletonGraph' in the simplicial complex. */ + template + void insert_graph(const OneSkeletonGraph& skel_graph); + + /** \brief Expands the simplicial complex containing only its one skeleton until a given maximal dimension. */ + void expansion(int max_dim); + +}; + +} // namespace rips_complex + +} // namespace Gudhi + +#endif // CONCEPT_RIPS_COMPLEX_SIMPLICIAL_COMPLEX_FOR_RIPS_H_ diff --git a/src/Rips_complex/doc/Intro_rips_complex.h b/src/Rips_complex/doc/Intro_rips_complex.h new file mode 100644 index 00000000..85108168 --- /dev/null +++ b/src/Rips_complex/doc/Intro_rips_complex.h @@ -0,0 +1,104 @@ +/* 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 & Vincent Rouvreau + * + * Copyright (C) 2015 INRIA + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ + +#ifndef DOC_RIPS_COMPLEX_INTRO_RIPS_COMPLEX_H_ +#define DOC_RIPS_COMPLEX_INTRO_RIPS_COMPLEX_H_ + +namespace Gudhi { + +namespace rips_complex { + +/** \defgroup rips_complex Rips complex + * + * \author Clément Maria & Vincent Rouvreau + * + * @{ + * + * \section definition Definition + * + * Rips_complex + * (Wikipedia) is a + * simplicial complex + * constructed from a one skeleton graph. + * + * The filtration value of each edge is computed from a user-given distance function. + * + * All edges that have a filtration value strictly greater than a given threshold value are not inserted into + * the complex. + * + * When creating a simplicial complex from this one skeleton graph, rips inserts the one skeleton graph into the data + * structure, and then expands the simplicial when required. + * + * \image html "rips_complex_representation.png" "Rips-complex one skeleton graph representation" + * + * On this example, as edges (4,5), (4,6) and (5,6) are in the complex, simplex (4,5,6) is added with the filtration + * value set with \f$max(filtration(4,5), filtration(4,6), filtration(5,6))\f$. + * And so on for simplex (0,1,2,3). + * + * \section ripspointsexample Example from points + * + * This example builds the one skeleton graph from the given points, threshold value, and distance function. + * Then it creates a `Simplex_tree` with it. + * + * Then, it is asked to display information about the simplicial complex. + * + * \include Rips_complex/example_rips_complex_from_points.cpp + * + * When launching: + * + * \code $> ./ripspoints + * \endcode + * + * the program output is: + * + * \include Rips_complex/rips_points_for_doc_12_2.txt + * + * \section offexample Example from OFF file + * + * This example builds the one skeleton graph from the given points in an OFF file, threshold value, and distance + * function. + * Then it creates a `Simplex_tree` with it. + * + * + * Then, it is asked to display information about the rips complex. + * + * \include Rips_complex/example_rips_complex_from_off_file.cpp + * + * When launching: + * + * \code $> ./ripsoffreader ../../data/points/alphacomplexdoc.off 12.0 3 + * \endcode + * + * the program output is: + * + * \include Rips_complex/rips_points_for_doc_12_3.txt + * + * \copyright GNU General Public License v3. + * \verbatim Contact: gudhi-users@lists.gforge.inria.fr \endverbatim + */ +/** @} */ // end defgroup rips_complex + +} // namespace rips_complex + +} // namespace Gudhi + +#endif // DOC_RIPS_COMPLEX_INTRO_RIPS_COMPLEX_H_ diff --git a/src/Rips_complex/doc/rips_complex_representation.ipe b/src/Rips_complex/doc/rips_complex_representation.ipe new file mode 100644 index 00000000..7f6028f4 --- /dev/null +++ b/src/Rips_complex/doc/rips_complex_representation.ipe @@ -0,0 +1,326 @@ + + + + + + + +0 0 m +-1 0.333 l +-1 -0.333 l +h + + + + +0 0 m +-1 0.333 l +-1 -0.333 l +h + + + + +0.6 0 0 0.6 0 0 e +0.4 0 0 0.4 0 0 e + + + + +0.6 0 0 0.6 0 0 e + + + + + +0.5 0 0 0.5 0 0 e + + +0.6 0 0 0.6 0 0 e +0.4 0 0 0.4 0 0 e + + + + + +-0.6 -0.6 m +0.6 -0.6 l +0.6 0.6 l +-0.6 0.6 l +h +-0.4 -0.4 m +0.4 -0.4 l +0.4 0.4 l +-0.4 0.4 l +h + + + + +-0.6 -0.6 m +0.6 -0.6 l +0.6 0.6 l +-0.6 0.6 l +h + + + + + +-0.5 -0.5 m +0.5 -0.5 l +0.5 0.5 l +-0.5 0.5 l +h + + +-0.6 -0.6 m +0.6 -0.6 l +0.6 0.6 l +-0.6 0.6 l +h +-0.4 -0.4 m +0.4 -0.4 l +0.4 0.4 l +-0.4 0.4 l +h + + + + + + +-0.43 -0.57 m +0.57 0.43 l +0.43 0.57 l +-0.57 -0.43 l +h + + +-0.43 0.57 m +0.57 -0.43 l +0.43 -0.57 l +-0.57 0.43 l +h + + + + + +0 0 m +-1 0.333 l +-1 -0.333 l +h + + + + +0 0 m +-1 0.333 l +-0.8 0 l +-1 -0.333 l +h + + + + +0 0 m +-1 0.333 l +-0.8 0 l +-1 -0.333 l +h + + + + +-1 0.333 m +0 0 l +-1 -0.333 l + + + + +0 0 m +-1 0.333 l +-1 -0.333 l +h +-1 0 m +-2 0.333 l +-2 -0.333 l +h + + + + +0 0 m +-1 0.333 l +-1 -0.333 l +h +-1 0 m +-2 0.333 l +-2 -0.333 l +h + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +109.771 601.912 m +159.595 601.797 l +140.058 541.915 l +h + + +79.8776 552.169 m +109.756 601.699 l +139.812 542.209 l +h + + +69.8453 682.419 m +159.925 712.208 l +90.12 732.039 l +h + +Rips complex +0 +1 +2 +3 +4 +5 +6 + +60 710 m +40 660 l + + +40 660 m +130 690 l + + +130 690 m +60 710 l + + +40 660 m +80 580 l + + +80 580 m +130 580 l +130 580 l + + +130 580 m +110 520 l + + +110 520 m +50 530 l +50 530 l +50 530 l + + +50 530 m +80 580 l + + +130 580 m +130 690 l + + + + + + + +150.038 609.9 m +179.929 549.727 l + + + + +158.7 593.269 m +81.4925 544.805 l + + +256.324 639.958 m +370.055 639.958 l + + +56.8567 0 0 56.8567 313.217 639.756 e + + + +Rips threshold + + diff --git a/src/Rips_complex/doc/rips_complex_representation.png b/src/Rips_complex/doc/rips_complex_representation.png new file mode 100644 index 00000000..e901d92e Binary files /dev/null and b/src/Rips_complex/doc/rips_complex_representation.png differ diff --git a/src/Rips_complex/doc/rips_one_skeleton.ipe b/src/Rips_complex/doc/rips_one_skeleton.ipe new file mode 100644 index 00000000..3a35970c --- /dev/null +++ b/src/Rips_complex/doc/rips_one_skeleton.ipe @@ -0,0 +1,326 @@ + + + + + + + +0 0 m +-1 0.333 l +-1 -0.333 l +h + + + + +0 0 m +-1 0.333 l +-1 -0.333 l +h + + + + +0.6 0 0 0.6 0 0 e +0.4 0 0 0.4 0 0 e + + + + +0.6 0 0 0.6 0 0 e + + + + + +0.5 0 0 0.5 0 0 e + + +0.6 0 0 0.6 0 0 e +0.4 0 0 0.4 0 0 e + + + + + +-0.6 -0.6 m +0.6 -0.6 l +0.6 0.6 l +-0.6 0.6 l +h +-0.4 -0.4 m +0.4 -0.4 l +0.4 0.4 l +-0.4 0.4 l +h + + + + +-0.6 -0.6 m +0.6 -0.6 l +0.6 0.6 l +-0.6 0.6 l +h + + + + + +-0.5 -0.5 m +0.5 -0.5 l +0.5 0.5 l +-0.5 0.5 l +h + + +-0.6 -0.6 m +0.6 -0.6 l +0.6 0.6 l +-0.6 0.6 l +h +-0.4 -0.4 m +0.4 -0.4 l +0.4 0.4 l +-0.4 0.4 l +h + + + + + + +-0.43 -0.57 m +0.57 0.43 l +0.43 0.57 l +-0.57 -0.43 l +h + + +-0.43 0.57 m +0.57 -0.43 l +0.43 -0.57 l +-0.57 0.43 l +h + + + + + +0 0 m +-1 0.333 l +-1 -0.333 l +h + + + + +0 0 m +-1 0.333 l +-0.8 0 l +-1 -0.333 l +h + + + + +0 0 m +-1 0.333 l +-0.8 0 l +-1 -0.333 l +h + + + + +-1 0.333 m +0 0 l +-1 -0.333 l + + + + +0 0 m +-1 0.333 l +-1 -0.333 l +h +-1 0 m +-2 0.333 l +-2 -0.333 l +h + + + + +0 0 m +-1 0.333 l +-1 -0.333 l +h +-1 0 m +-2 0.333 l +-2 -0.333 l +h + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +109.771 601.912 m +159.595 601.797 l +140.058 541.915 l +h + + +79.8776 552.169 m +109.756 601.699 l +139.812 542.209 l +h + + +69.8453 682.419 m +159.925 712.208 l +90.12 732.039 l +h + +One skeleton graph +0 +1 +2 +3 +4 +5 +6 + +60 710 m +40 660 l + + +40 660 m +130 690 l + + +130 690 m +60 710 l + + +40 660 m +80 580 l + + +80 580 m +130 580 l +130 580 l + + +130 580 m +110 520 l + + +110 520 m +50 530 l +50 530 l +50 530 l + + +50 530 m +80 580 l + + +130 580 m +130 690 l + + + + + + + +150.038 609.9 m +179.929 549.727 l + + + + +158.7 593.269 m +81.4925 544.805 l + + +256.324 639.958 m +370.055 639.958 l + + +56.8567 0 0 56.8567 313.217 639.756 e + + + +Rips threshold + + diff --git a/src/Rips_complex/doc/rips_one_skeleton.png b/src/Rips_complex/doc/rips_one_skeleton.png new file mode 100644 index 00000000..1028770e Binary files /dev/null and b/src/Rips_complex/doc/rips_one_skeleton.png differ diff --git a/src/Rips_complex/example/CMakeLists.txt b/src/Rips_complex/example/CMakeLists.txt new file mode 100644 index 00000000..6d0deecf --- /dev/null +++ b/src/Rips_complex/example/CMakeLists.txt @@ -0,0 +1,8 @@ +cmake_minimum_required(VERSION 2.6) +project(Rips_complex_examples) + +add_executable ( ripsoffreader example_rips_complex_from_off_file.cpp ) +target_link_libraries(ripsoffreader ${Boost_SYSTEM_LIBRARY}) +if (TBB_FOUND) + target_link_libraries(ripsoffreader ${TBB_LIBRARIES}) +endif() diff --git a/src/Rips_complex/example/example_rips_complex_from_off_file.cpp b/src/Rips_complex/example/example_rips_complex_from_off_file.cpp new file mode 100644 index 00000000..82baa68e --- /dev/null +++ b/src/Rips_complex/example/example_rips_complex_from_off_file.cpp @@ -0,0 +1,71 @@ +#include +// to construct Rips_complex from a OFF file of points +#include +// to construct a simplex_tree from rips complex +#include +#include + +#include +#include + +void usage(int nbArgs, char * const progName) { + std::cerr << "Error: Number of arguments (" << nbArgs << ") is not correct\n"; + std::cerr << "Usage: " << progName << " filename.off threshold dim_max [ouput_file.txt]\n"; + std::cerr << " i.e.: " << progName << " ../../data/points/alphacomplexdoc.off 60.0\n"; + exit(-1); // ----- >> +} + +int main(int argc, char **argv) { + if ((argc != 4) && (argc != 5)) usage(argc, (argv[0] - 1)); + + std::string off_file_name(argv[1]); + double threshold = atof(argv[2]); + int dim_max = atoi(argv[3]); + + // Type definitions + using Point = std::vector; + using Simplex_tree = Gudhi::Simplex_tree<>; + using Rips_complex = Gudhi::rips_complex::Rips_complex; + + // ---------------------------------------------------------------------------- + // Init of a rips complex from an OFF file + // ---------------------------------------------------------------------------- + Gudhi::Points_off_reader off_reader(off_file_name); + Rips_complex rips_complex_from_file(off_reader.get_point_cloud(), threshold, + euclidean_distance); + + std::streambuf* streambufffer; + std::ofstream ouput_file_stream; + + if (argc == 5) { + ouput_file_stream.open(std::string(argv[4])); + streambufffer = ouput_file_stream.rdbuf(); + } else { + streambufffer = std::cout.rdbuf(); + } + + Simplex_tree simplex; + if (rips_complex_from_file.create_complex(simplex, dim_max)) { + std::ostream output_stream(streambufffer); + + // ---------------------------------------------------------------------------- + // Display information about the rips complex + // ---------------------------------------------------------------------------- + output_stream << "Rips complex is of dimension " << simplex.dimension() << + " - " << simplex.num_simplices() << " simplices - " << + simplex.num_vertices() << " vertices." << std::endl; + + output_stream << "Iterator on rips complex simplices in the filtration order, with [filtration value]:" << + std::endl; + for (auto f_simplex : simplex.filtration_simplex_range()) { + output_stream << " ( "; + for (auto vertex : simplex.simplex_vertex_range(f_simplex)) { + output_stream << vertex << " "; + } + output_stream << ") -> " << "[" << simplex.filtration(f_simplex) << "] "; + output_stream << std::endl; + } + } + ouput_file_stream.close(); + return 0; +} diff --git a/src/Rips_complex/include/gudhi/Rips_complex.h b/src/Rips_complex/include/gudhi/Rips_complex.h new file mode 100644 index 00000000..10f674e5 --- /dev/null +++ b/src/Rips_complex/include/gudhi/Rips_complex.h @@ -0,0 +1,154 @@ +/* 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 & Vincent Rouvreau + * + * Copyright (C) 2016 INRIA + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ + +#ifndef RIPS_COMPLEX_H_ +#define RIPS_COMPLEX_H_ + +#include +#include + +#include + +#include +#include +#include +#include +#include // for numeric_limits +#include // for pair<> + + +namespace Gudhi { + +namespace rips_complex { + +/** + * \class Rips_complex + * \brief Rips complex data structure. + * + * \ingroup rips_complex + * + * \details + * The data structure is a 1-skeleton graph constructed from a point cloud, containing edges when the edge length is + * less or equal to a given threshold. Edge length is computed from a user given function. + * + * The complex is a template class requiring a Filtration_value type. + * + * \remark When Alpha_complex is constructed with an infinite value of alpha, the complex is a Delaunay complex. + * + * \tparam Filtration_value must meet `SimplicialComplexForRips` concept. + */ +template +class Rips_complex { + private: + typedef typename boost::adjacency_list < boost::vecS, boost::vecS, boost::undirectedS + , boost::property < vertex_filtration_t, Filtration_value > + , boost::property < edge_filtration_t, Filtration_value >> Graph_t; + + typedef int Vertex_handle; + + public: + /** \brief Rips_complex constructor from a list of points. + * + * @param[in] points Range of points. + * @param[in] threshold rips value. + * @param[in] distance distance function that returns a Filtration_value from 2 given points. + * + * The type InputPointRange must be a range for which std::begin and std::end return input iterators on a point. + */ + template + Rips_complex(const InputPointRange& points, Filtration_value threshold, + Filtration_value distance(const Point_d& p1,const Point_d& p2)) { + std::vector< std::pair< Vertex_handle, Vertex_handle > > edges; + std::vector< Filtration_value > edges_fil; + std::map< Vertex_handle, Filtration_value > vertices; + + // Compute the proximity graph of the points. + // If points contains n elements, the proximity graph is the graph with n vertices, and an edge [u,v] iff the + // distance function between points u and v is smaller than threshold. + // -------------------------------------------------------------------------------------------- + // Creates the vector of edges and its filtration values (returned by distance function) + Vertex_handle idx_u, idx_v; + Filtration_value fil; + idx_u = 0; + for (auto it_u = std::begin(points); it_u != std::end(points); ++it_u) { + idx_v = idx_u + 1; + for (auto it_v = it_u + 1; it_v != std::end(points); ++it_v, ++idx_v) { + fil = distance(*it_u, *it_v); + if (fil <= threshold) { + edges.emplace_back(idx_u, idx_v); + edges_fil.push_back(fil); + } + } + ++idx_u; + } + + // -------------------------------------------------------------------------------------------- + // Creates the proximity graph from edges and sets the property with the filtration value. + // Number of points is labeled from 0 to idx_u-1 + rips_skeleton_graph_ = Graph_t(edges.begin() , edges.end() , edges_fil.begin() , idx_u); + + auto vertex_prop = boost::get(vertex_filtration_t(), rips_skeleton_graph_); + + using vertex_iterator = typename boost::graph_traits::vertex_iterator; + vertex_iterator vi, vi_end; + for (std::tie(vi, vi_end) = boost::vertices(rips_skeleton_graph_); + vi != vi_end; ++vi) { + boost::put(vertex_prop, *vi, 0.); + } + + } + + /** \brief Initializes the simplicial complex from the 1-skeleton graph and expands it until a given maximal + * dimension. + * + * \tparam SimplicialComplexForRips must meet `SimplicialComplexForRips` concept. + * + * @param[in] complex SimplicialComplexForRips to be created. + * @param[in] dim_max graph expansion for rips until this given maximal dimension. + * + * @return true if creation succeeds, false otherwise. + * + */ + template + bool create_complex(SimplicialComplexForRips& complex, int dim_max) { + if (complex.num_vertices() > 0) { + std::cerr << "Rips_complex create_complex - complex is not empty\n"; + return false; // ----- >> + } + + // insert the proximity graph in the simplicial complex + complex.insert_graph(rips_skeleton_graph_); + // expand the graph until dimension dim_max + complex.expansion(dim_max); + + // -------------------------------------------------------------------------------------------- + return true; + } + private: + Graph_t rips_skeleton_graph_; +}; + +} // namespace rips_complex + +} // namespace Gudhi + +#endif // RIPS_COMPLEX_H_ diff --git a/src/Simplex_tree/include/gudhi/Simplex_tree.h b/src/Simplex_tree/include/gudhi/Simplex_tree.h index 63e3f0e5..206b2fba 100644 --- a/src/Simplex_tree/include/gudhi/Simplex_tree.h +++ b/src/Simplex_tree/include/gudhi/Simplex_tree.h @@ -28,7 +28,7 @@ #include #include -#include +//#include #include #include diff --git a/src/common/include/gudhi/distance_functions.h b/src/common/include/gudhi/distance_functions.h index cd518581..b2726ba8 100644 --- a/src/common/include/gudhi/distance_functions.h +++ b/src/common/include/gudhi/distance_functions.h @@ -29,7 +29,7 @@ * by a range of coordinates. The points are assumed to have * the same dimension. */ template< typename Point > -double euclidean_distance(Point &p1, Point &p2) { +double euclidean_distance(const Point &p1,const Point &p2) { double dist = 0.; auto it1 = p1.begin(); auto it2 = p2.begin(); diff --git a/src/common/include/gudhi/graph_simplicial_complex.h b/src/common/include/gudhi/graph_simplicial_complex.h index 042ef516..773889d9 100644 --- a/src/common/include/gudhi/graph_simplicial_complex.h +++ b/src/common/include/gudhi/graph_simplicial_complex.h @@ -39,14 +39,14 @@ struct vertex_filtration_t { typedef boost::vertex_property_tag kind; }; -typedef int Vertex_handle; +/*typedef int Vertex_handle; typedef double Filtration_value; typedef boost::adjacency_list < boost::vecS, boost::vecS, boost::undirectedS , boost::property < vertex_filtration_t, Filtration_value > , boost::property < edge_filtration_t, Filtration_value > > Graph_t; typedef std::pair< Vertex_handle, Vertex_handle > Edge_t; - +*/ /** \brief Output the proximity graph of the points. * * If points contains n elements, the proximity graph is the graph @@ -56,7 +56,7 @@ typedef std::pair< Vertex_handle, Vertex_handle > Edge_t; * The type PointCloud furnishes .begin() and .end() methods, that return * iterators with value_type Point. */ -template< typename PointCloud +/*template< typename PointCloud , typename Point > Graph_t compute_proximity_graph(PointCloud &points , Filtration_value threshold @@ -94,6 +94,6 @@ Graph_t compute_proximity_graph(PointCloud &points } return skel_graph; -} +}*/ #endif // GRAPH_SIMPLICIAL_COMPLEX_H_ -- cgit v1.2.3 From 1c1c3073e0612ee401a941008a6b13438b3429e2 Mon Sep 17 00:00:00 2001 From: vrouvrea Date: Thu, 29 Sep 2016 12:34:19 +0000 Subject: Fix doxygen Add tests git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/rips_complex_module@1589 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: 3b03e969f43f4f41c7bca1b2de96a71688475742 --- CMakeLists.txt | 1 + src/CMakeLists.txt | 1 + src/Doxyfile | 3 +- src/Rips_complex/doc/Intro_rips_complex.h | 17 +- src/Rips_complex/example/CMakeLists.txt | 21 +- .../example_one_skeleton_rips_from_points.cpp | 64 ++++++ .../example/example_rips_complex_from_points.cpp | 69 ------ .../example/full_skeleton_rips_points_for_doc.txt | 26 +++ .../example/one_skeleton_rips_points_for_doc.txt | 20 ++ .../example/rips_points_for_doc_12_1.txt | 20 -- .../example/rips_points_for_doc_12_3.txt | 26 --- src/Rips_complex/test/CMakeLists.txt | 24 +++ src/Rips_complex/test/README | 12 ++ src/Rips_complex/test/test_rips_complex.cpp | 231 +++++++++++++++++++++ src/cmake/modules/GUDHI_user_version_target.txt | 2 +- src/common/doc/main_page.h | 23 ++ 16 files changed, 431 insertions(+), 129 deletions(-) create mode 100644 src/Rips_complex/example/example_one_skeleton_rips_from_points.cpp delete mode 100644 src/Rips_complex/example/example_rips_complex_from_points.cpp create mode 100644 src/Rips_complex/example/full_skeleton_rips_points_for_doc.txt create mode 100644 src/Rips_complex/example/one_skeleton_rips_points_for_doc.txt delete mode 100644 src/Rips_complex/example/rips_points_for_doc_12_1.txt delete mode 100644 src/Rips_complex/example/rips_points_for_doc_12_3.txt create mode 100644 src/Rips_complex/test/CMakeLists.txt create mode 100644 src/Rips_complex/test/README create mode 100644 src/Rips_complex/test/test_rips_complex.cpp (limited to 'CMakeLists.txt') diff --git a/CMakeLists.txt b/CMakeLists.txt index 89e130f4..972df204 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -144,6 +144,7 @@ else() add_subdirectory(src/Alpha_complex/example) add_subdirectory(src/Alpha_complex/test) add_subdirectory(src/Rips_complex/example) + add_subdirectory(src/Rips_complex/test) # data points generator add_subdirectory(data/points/generator) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index c7744c49..7973ae39 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -108,6 +108,7 @@ else() add_subdirectory(example/Bitmap_cubical_complex) add_subdirectory(example/Witness_complex) add_subdirectory(example/Alpha_complex) + add_subdirectory(example/Rips_complex) # data points generator add_subdirectory(data/points/generator) diff --git a/src/Doxyfile b/src/Doxyfile index dd9a33fb..40d07049 100644 --- a/src/Doxyfile +++ b/src/Doxyfile @@ -845,7 +845,8 @@ IMAGE_PATH = doc/Skeleton_blocker/ \ doc/Simplex_tree/ \ doc/Persistent_cohomology/ \ doc/Witness_complex/ \ - doc/Bitmap_cubical_complex/ + doc/Bitmap_cubical_complex/ \ + doc/Rips_complex/ # The INPUT_FILTER tag can be used to specify a program that doxygen should # invoke to filter for each input file. Doxygen will invoke the filter program diff --git a/src/Rips_complex/doc/Intro_rips_complex.h b/src/Rips_complex/doc/Intro_rips_complex.h index 85108168..bd9f5fad 100644 --- a/src/Rips_complex/doc/Intro_rips_complex.h +++ b/src/Rips_complex/doc/Intro_rips_complex.h @@ -29,11 +29,11 @@ namespace rips_complex { /** \defgroup rips_complex Rips complex * - * \author Clément Maria & Vincent Rouvreau + * \author Clément Maria and Vincent Rouvreau * * @{ * - * \section definition Definition + * \section ripsdefinition Definition * * Rips_complex * (Wikipedia) is a @@ -61,18 +61,19 @@ namespace rips_complex { * * Then, it is asked to display information about the simplicial complex. * - * \include Rips_complex/example_rips_complex_from_points.cpp + * \include Rips_complex/example_one_skeleton_rips_from_points.cpp * - * When launching: + * When launching (rips maximal distance between 2 points is 12.0, is expanded until dimension 1 - one skeleton graph + * with other words): * - * \code $> ./ripspoints + * \code $> ./oneskeletonripspoints 12.0 * \endcode * * the program output is: * - * \include Rips_complex/rips_points_for_doc_12_2.txt + * \include Rips_complex/one_skeleton_rips_points_for_doc.txt * - * \section offexample Example from OFF file + * \section ripsoffexample Example from OFF file * * This example builds the one skeleton graph from the given points in an OFF file, threshold value, and distance * function. @@ -90,7 +91,7 @@ namespace rips_complex { * * the program output is: * - * \include Rips_complex/rips_points_for_doc_12_3.txt + * \include Rips_complex/full_skeleton_rips_points_for_doc.txt * * \copyright GNU General Public License v3. * \verbatim Contact: gudhi-users@lists.gforge.inria.fr \endverbatim diff --git a/src/Rips_complex/example/CMakeLists.txt b/src/Rips_complex/example/CMakeLists.txt index f8dcf94d..3962b15e 100644 --- a/src/Rips_complex/example/CMakeLists.txt +++ b/src/Rips_complex/example/CMakeLists.txt @@ -3,12 +3,25 @@ project(Rips_complex_examples) add_executable ( ripsoffreader example_rips_complex_from_off_file.cpp ) target_link_libraries(ripsoffreader ${Boost_SYSTEM_LIBRARY}) + +add_executable ( oneskeletonripspoints example_one_skeleton_rips_from_points.cpp ) +target_link_libraries(oneskeletonripspoints ${Boost_SYSTEM_LIBRARY}) + if (TBB_FOUND) target_link_libraries(ripsoffreader ${TBB_LIBRARIES}) + target_link_libraries(oneskeletonripspoints ${TBB_LIBRARIES}) endif() -add_executable ( ripspoints example_rips_complex_from_points.cpp ) -target_link_libraries(ripspoints ${Boost_SYSTEM_LIBRARY}) -if (TBB_FOUND) - target_link_libraries(ripspoints ${TBB_LIBRARIES}) +add_test(oneskeletonripspoints ${CMAKE_CURRENT_BINARY_DIR}/oneskeletonripspoints 12.0) +# Do not forget to copy test files in current binary dir +file(COPY "${CMAKE_SOURCE_DIR}/data/points/alphacomplexdoc.off" DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/) +add_test(ripsoffreader_doc_12_1 ${CMAKE_CURRENT_BINARY_DIR}/ripsoffreader alphacomplexdoc.off 12.0 1 ${CMAKE_CURRENT_BINARY_DIR}/ripsoffreader_result_12_1.txt) +add_test(ripsoffreader_doc_12_3 ${CMAKE_CURRENT_BINARY_DIR}/ripsoffreader alphacomplexdoc.off 12.0 3 ${CMAKE_CURRENT_BINARY_DIR}/ripsoffreader_result_12_3.txt) +if (DIFF_PATH) + # Do not forget to copy test results files in current binary dir + file(COPY "one_skeleton_rips_points_for_doc.txt" DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/) + file(COPY "full_skeleton_rips_points_for_doc.txt" DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/) + + add_test(ripsoffreader_doc_12_1_diff_files ${DIFF_PATH} ${CMAKE_CURRENT_BINARY_DIR}/ripsoffreader_result_12_1.txt ${CMAKE_CURRENT_BINARY_DIR}/one_skeleton_rips_points_for_doc.txt) + add_test(ripsoffreader_doc_12_3_diff_files ${DIFF_PATH} ${CMAKE_CURRENT_BINARY_DIR}/ripsoffreader_result_12_3.txt ${CMAKE_CURRENT_BINARY_DIR}/full_skeleton_rips_points_for_doc.txt) endif() diff --git a/src/Rips_complex/example/example_one_skeleton_rips_from_points.cpp b/src/Rips_complex/example/example_one_skeleton_rips_from_points.cpp new file mode 100644 index 00000000..db41b425 --- /dev/null +++ b/src/Rips_complex/example/example_one_skeleton_rips_from_points.cpp @@ -0,0 +1,64 @@ +#include +// to construct Rips_complex from a OFF file of points +#include +// to construct a simplex_tree from rips complex +#include +#include + +#include +#include +#include // for std::numeric_limits + +void usage(int nbArgs, char * const progName) { + std::cerr << "Error: Number of arguments (" << nbArgs << ") is not correct\n"; + std::cerr << "Usage: " << progName << " threshold\n"; + std::cerr << " i.e.: " << progName << " 12.0\n"; + exit(-1); // ----- >> +} + +int main(int argc, char **argv) { + if (argc != 2) usage(argc, argv[0]); + + double threshold = atof(argv[1]); + + // Type definitions + using Point = std::vector; + using Simplex_tree = Gudhi::Simplex_tree<>; + using Rips_complex = Gudhi::rips_complex::Rips_complex; + + std::vector points; + points.push_back({1.0, 1.0}); + points.push_back({7.0, 0.0}); + points.push_back({4.0, 6.0}); + points.push_back({9.0, 6.0}); + points.push_back({0.0, 14.0}); + points.push_back({2.0, 19.0}); + points.push_back({9.0, 17.0}); + + // ---------------------------------------------------------------------------- + // Init of a rips complex from points + // ---------------------------------------------------------------------------- + Rips_complex rips_complex_from_file(points, threshold, euclidean_distance); + + Simplex_tree simplex; + if (rips_complex_from_file.create_complex(simplex, 1)) { + // ---------------------------------------------------------------------------- + // Display information about the one skeleton rips complex + // ---------------------------------------------------------------------------- + std::cout << "Rips complex is of dimension " << simplex.dimension() << + " - " << simplex.num_simplices() << " simplices - " << + simplex.num_vertices() << " vertices." << std::endl; + + std::cout << "Iterator on rips complex simplices in the filtration order, with [filtration value]:" << + std::endl; + for (auto f_simplex : simplex.filtration_simplex_range()) { + std::cout << " ( "; + for (auto vertex : simplex.simplex_vertex_range(f_simplex)) { + std::cout << vertex << " "; + } + std::cout << ") -> " << "[" << simplex.filtration(f_simplex) << "] "; + std::cout << std::endl; + } + } + return 0; +} diff --git a/src/Rips_complex/example/example_rips_complex_from_points.cpp b/src/Rips_complex/example/example_rips_complex_from_points.cpp deleted file mode 100644 index c241d137..00000000 --- a/src/Rips_complex/example/example_rips_complex_from_points.cpp +++ /dev/null @@ -1,69 +0,0 @@ -#include -// to construct Rips_complex from a OFF file of points -#include -// to construct a simplex_tree from rips complex -#include -#include - -#include -#include -#include // for std::numeric_limits - -void usage(int nbArgs, char * const progName) { - std::cerr << "Error: Number of arguments (" << nbArgs << ") is not correct\n"; - std::cerr << "Usage: " << progName << " dim_max [threshold]\n"; - std::cerr << " i.e.: " << progName << " 3 12.0\n"; - exit(-1); // ----- >> -} - -int main(int argc, char **argv) { - if ((argc != 2) && (argc != 3)) usage(argc, argv[0]); - - double threshold = std::numeric_limits::infinity(); - int dim_max = atoi(argv[1]); - - if (argc == 3) { - threshold = atof(argv[2]); - } - - // Type definitions - using Point = std::vector; - using Simplex_tree = Gudhi::Simplex_tree<>; - using Rips_complex = Gudhi::rips_complex::Rips_complex; - - std::vector points; - points.push_back({1.0, 1.0}); - points.push_back({7.0, 0.0}); - points.push_back({4.0, 6.0}); - points.push_back({9.0, 6.0}); - points.push_back({0.0, 14.0}); - points.push_back({2.0, 19.0}); - points.push_back({9.0, 17.0}); - - // ---------------------------------------------------------------------------- - // Init of a rips complex from an OFF file - // ---------------------------------------------------------------------------- - Rips_complex rips_complex_from_file(points, threshold, euclidean_distance); - - Simplex_tree simplex; - if (rips_complex_from_file.create_complex(simplex, dim_max)) { - // ---------------------------------------------------------------------------- - // Display information about the rips complex - // ---------------------------------------------------------------------------- - std::cout << "Rips complex is of dimension " << simplex.dimension() << - " - " << simplex.num_simplices() << " simplices - " << - simplex.num_vertices() << " vertices." << std::endl; - - std::cout << "Iterator on rips complex simplices in the filtration order, with [filtration value]:" << - std::endl; - for (auto f_simplex : simplex.filtration_simplex_range()) { - std::cout << " ( "; - for (auto vertex : simplex.simplex_vertex_range(f_simplex)) { - std::cout << vertex << " "; - } - std::cout << ") -> " << "[" << simplex.filtration(f_simplex) << "] "; - std::cout << std::endl; - } - } - return 0; -} diff --git a/src/Rips_complex/example/full_skeleton_rips_points_for_doc.txt b/src/Rips_complex/example/full_skeleton_rips_points_for_doc.txt new file mode 100644 index 00000000..319931e0 --- /dev/null +++ b/src/Rips_complex/example/full_skeleton_rips_points_for_doc.txt @@ -0,0 +1,26 @@ +Rips complex is of dimension 3 - 24 simplices - 7 vertices. +Iterator on rips complex simplices in the filtration order, with [filtration value]: + ( 0 ) -> [0] + ( 1 ) -> [0] + ( 2 ) -> [0] + ( 3 ) -> [0] + ( 4 ) -> [0] + ( 5 ) -> [0] + ( 6 ) -> [0] + ( 3 2 ) -> [5] + ( 5 4 ) -> [5.38516] + ( 2 0 ) -> [5.83095] + ( 1 0 ) -> [6.08276] + ( 3 1 ) -> [6.32456] + ( 2 1 ) -> [6.7082] + ( 2 1 0 ) -> [6.7082] + ( 3 2 1 ) -> [6.7082] + ( 6 5 ) -> [7.28011] + ( 4 2 ) -> [8.94427] + ( 3 0 ) -> [9.43398] + ( 3 1 0 ) -> [9.43398] + ( 3 2 0 ) -> [9.43398] + ( 3 2 1 0 ) -> [9.43398] + ( 6 4 ) -> [9.48683] + ( 6 5 4 ) -> [9.48683] + ( 6 3 ) -> [11] diff --git a/src/Rips_complex/example/one_skeleton_rips_points_for_doc.txt b/src/Rips_complex/example/one_skeleton_rips_points_for_doc.txt new file mode 100644 index 00000000..b0e25cc5 --- /dev/null +++ b/src/Rips_complex/example/one_skeleton_rips_points_for_doc.txt @@ -0,0 +1,20 @@ +Rips complex is of dimension 1 - 18 simplices - 7 vertices. +Iterator on rips complex simplices in the filtration order, with [filtration value]: + ( 0 ) -> [0] + ( 1 ) -> [0] + ( 2 ) -> [0] + ( 3 ) -> [0] + ( 4 ) -> [0] + ( 5 ) -> [0] + ( 6 ) -> [0] + ( 3 2 ) -> [5] + ( 5 4 ) -> [5.38516] + ( 2 0 ) -> [5.83095] + ( 1 0 ) -> [6.08276] + ( 3 1 ) -> [6.32456] + ( 2 1 ) -> [6.7082] + ( 6 5 ) -> [7.28011] + ( 4 2 ) -> [8.94427] + ( 3 0 ) -> [9.43398] + ( 6 4 ) -> [9.48683] + ( 6 3 ) -> [11] diff --git a/src/Rips_complex/example/rips_points_for_doc_12_1.txt b/src/Rips_complex/example/rips_points_for_doc_12_1.txt deleted file mode 100644 index b0e25cc5..00000000 --- a/src/Rips_complex/example/rips_points_for_doc_12_1.txt +++ /dev/null @@ -1,20 +0,0 @@ -Rips complex is of dimension 1 - 18 simplices - 7 vertices. -Iterator on rips complex simplices in the filtration order, with [filtration value]: - ( 0 ) -> [0] - ( 1 ) -> [0] - ( 2 ) -> [0] - ( 3 ) -> [0] - ( 4 ) -> [0] - ( 5 ) -> [0] - ( 6 ) -> [0] - ( 3 2 ) -> [5] - ( 5 4 ) -> [5.38516] - ( 2 0 ) -> [5.83095] - ( 1 0 ) -> [6.08276] - ( 3 1 ) -> [6.32456] - ( 2 1 ) -> [6.7082] - ( 6 5 ) -> [7.28011] - ( 4 2 ) -> [8.94427] - ( 3 0 ) -> [9.43398] - ( 6 4 ) -> [9.48683] - ( 6 3 ) -> [11] diff --git a/src/Rips_complex/example/rips_points_for_doc_12_3.txt b/src/Rips_complex/example/rips_points_for_doc_12_3.txt deleted file mode 100644 index 319931e0..00000000 --- a/src/Rips_complex/example/rips_points_for_doc_12_3.txt +++ /dev/null @@ -1,26 +0,0 @@ -Rips complex is of dimension 3 - 24 simplices - 7 vertices. -Iterator on rips complex simplices in the filtration order, with [filtration value]: - ( 0 ) -> [0] - ( 1 ) -> [0] - ( 2 ) -> [0] - ( 3 ) -> [0] - ( 4 ) -> [0] - ( 5 ) -> [0] - ( 6 ) -> [0] - ( 3 2 ) -> [5] - ( 5 4 ) -> [5.38516] - ( 2 0 ) -> [5.83095] - ( 1 0 ) -> [6.08276] - ( 3 1 ) -> [6.32456] - ( 2 1 ) -> [6.7082] - ( 2 1 0 ) -> [6.7082] - ( 3 2 1 ) -> [6.7082] - ( 6 5 ) -> [7.28011] - ( 4 2 ) -> [8.94427] - ( 3 0 ) -> [9.43398] - ( 3 1 0 ) -> [9.43398] - ( 3 2 0 ) -> [9.43398] - ( 3 2 1 0 ) -> [9.43398] - ( 6 4 ) -> [9.48683] - ( 6 5 4 ) -> [9.48683] - ( 6 3 ) -> [11] diff --git a/src/Rips_complex/test/CMakeLists.txt b/src/Rips_complex/test/CMakeLists.txt new file mode 100644 index 00000000..3972e9fc --- /dev/null +++ b/src/Rips_complex/test/CMakeLists.txt @@ -0,0 +1,24 @@ +cmake_minimum_required(VERSION 2.6) +project(Rips_complex_tests) + +if (GCOVR_PATH) + # for gcovr to make coverage reports - Corbera Jenkins plugin + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -fprofile-arcs -ftest-coverage") +endif() +if (GPROF_PATH) + # for gprof to make coverage reports - Jenkins + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -pg") +endif() + +add_executable ( rips_complex_UT test_rips_complex.cpp ) +target_link_libraries(rips_complex_UT ${Boost_SYSTEM_LIBRARY} ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY}) +if (TBB_FOUND) + target_link_libraries(rips_complex_UT ${TBB_LIBRARIES}) +endif() + +# Do not forget to copy test files in current binary dir +file(COPY "${CMAKE_SOURCE_DIR}/data/points/alphacomplexdoc.off" DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/) + +add_test(rips_complex_UT ${CMAKE_CURRENT_BINARY_DIR}/rips_complex_UT + # XML format for Jenkins xUnit plugin + --log_format=XML --log_sink=${CMAKE_SOURCE_DIR}/rips_complex_UT.xml --log_level=test_suite --report_level=no) diff --git a/src/Rips_complex/test/README b/src/Rips_complex/test/README new file mode 100644 index 00000000..28236b52 --- /dev/null +++ b/src/Rips_complex/test/README @@ -0,0 +1,12 @@ +To compile: +*********** + +cmake . +make + +To launch with details: +*********************** + +./rips_complex_UT --report_level=detailed --log_level=all + + ==> echo $? returns 0 in case of success (non-zero otherwise) diff --git a/src/Rips_complex/test/test_rips_complex.cpp b/src/Rips_complex/test/test_rips_complex.cpp new file mode 100644 index 00000000..5f129160 --- /dev/null +++ b/src/Rips_complex/test/test_rips_complex.cpp @@ -0,0 +1,231 @@ +/* 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): Vincent Rouvreau + * + * Copyright (C) 2016 INRIA Saclay (France) + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ + +#define BOOST_TEST_DYN_LINK +#define BOOST_TEST_MODULE "rips_complex" +#include + +#include // float comparison +#include +#include +#include +#include // std::max + +#include +// to construct Rips_complex from a OFF file of points +#include +// to construct a simplex_tree from rips complex +#include +#include + +// Type definitions +using Point = std::vector; +using Simplex_tree = Gudhi::Simplex_tree<>; +using Rips_complex = Gudhi::rips_complex::Rips_complex; + +bool are_almost_the_same(float a, float b) { + return std::fabs(a - b) < std::numeric_limits::epsilon(); +} + +BOOST_AUTO_TEST_CASE(RIPS_DOC_OFF_file) { + // ---------------------------------------------------------------------------- + // + // Init of a rips complex from a OFF file + // + // ---------------------------------------------------------------------------- + std::string off_file_name("alphacomplexdoc.off"); + double rips_threshold = 12.0; + std::cout << "========== OFF FILE NAME = " << off_file_name << " - rips threshold=" << + rips_threshold << "==========" << std::endl; + + Gudhi::Points_off_reader off_reader(off_file_name); + Rips_complex rips_complex_from_file(off_reader.get_point_cloud(), rips_threshold, euclidean_distance); + + const int DIMENSION_1 = 1; + Simplex_tree st; + BOOST_CHECK(rips_complex_from_file.create_complex(st, DIMENSION_1)); + std::cout << "st.dimension()=" << st.dimension() << std::endl; + BOOST_CHECK(st.dimension() == DIMENSION_1); + + const int NUMBER_OF_VERTICES = 7; + std::cout << "st.num_vertices()=" << st.num_vertices() << std::endl; + BOOST_CHECK(st.num_vertices() == NUMBER_OF_VERTICES); + + std::cout << "st.num_simplices()=" << st.num_simplices() << std::endl; + BOOST_CHECK(st.num_simplices() == 18); + + // Check filtration values of vertices is 0.0 + for (auto f_simplex : st.skeleton_simplex_range(0)) { + BOOST_CHECK(st.filtration(f_simplex) == 0.0); + } + + // Check filtration values of edges + for (auto f_simplex : st.skeleton_simplex_range(DIMENSION_1)) { + if (DIMENSION_1 == st.dimension(f_simplex)) { + std::vector vp; + std::cout << "vertex = ("; + for (auto vertex : st.simplex_vertex_range(f_simplex)) { + std::cout << vertex << ","; + vp.push_back(off_reader.get_point_cloud().at(vertex)); + } + std::cout << ") - distance =" << euclidean_distance(vp.at(0), vp.at(1)) << + " - filtration =" << st.filtration(f_simplex) << std::endl; + BOOST_CHECK(vp.size() == 2); + BOOST_CHECK(are_almost_the_same(st.filtration(f_simplex), euclidean_distance(vp.at(0), vp.at(1)))); + } + } + + const int DIMENSION_2 = 2; + Simplex_tree st2; + BOOST_CHECK(rips_complex_from_file.create_complex(st2, DIMENSION_2)); + std::cout << "st2.dimension()=" << st2.dimension() << std::endl; + BOOST_CHECK(st2.dimension() == DIMENSION_2); + + std::cout << "st2.num_vertices()=" << st2.num_vertices() << std::endl; + BOOST_CHECK(st2.num_vertices() == NUMBER_OF_VERTICES); + + std::cout << "st2.num_simplices()=" << st2.num_simplices() << std::endl; + BOOST_CHECK(st2.num_simplices() == 23); + + Simplex_tree::Filtration_value f01 = st2.filtration(st2.find({0, 1})); + Simplex_tree::Filtration_value f02 = st2.filtration(st2.find({0, 2})); + Simplex_tree::Filtration_value f12 = st2.filtration(st2.find({1, 2})); + Simplex_tree::Filtration_value f012 = st2.filtration(st2.find({0, 1, 2})); + std::cout << "f012= " << f012 << " | f01= " << f01 << " - f02= " << f02 << " - f12= " << f12 << std::endl; + BOOST_CHECK(are_almost_the_same(f012, std::max(f01, std::max(f02,f12)))); + + Simplex_tree::Filtration_value f45 = st2.filtration(st2.find({4, 5})); + Simplex_tree::Filtration_value f56 = st2.filtration(st2.find({5, 6})); + Simplex_tree::Filtration_value f46 = st2.filtration(st2.find({4, 6})); + Simplex_tree::Filtration_value f456 = st2.filtration(st2.find({4, 5, 6})); + std::cout << "f456= " << f456 << " | f45= " << f45 << " - f56= " << f56 << " - f46= " << f46 << std::endl; + BOOST_CHECK(are_almost_the_same(f456, std::max(f45, std::max(f56,f46)))); + + const int DIMENSION_3 = 3; + Simplex_tree st3; + BOOST_CHECK(rips_complex_from_file.create_complex(st3, DIMENSION_3)); + std::cout << "st3.dimension()=" << st3.dimension() << std::endl; + BOOST_CHECK(st3.dimension() == DIMENSION_3); + + std::cout << "st3.num_vertices()=" << st3.num_vertices() << std::endl; + BOOST_CHECK(st3.num_vertices() == NUMBER_OF_VERTICES); + + std::cout << "st3.num_simplices()=" << st3.num_simplices() << std::endl; + BOOST_CHECK(st3.num_simplices() == 24); + + Simplex_tree::Filtration_value f123 = st3.filtration(st3.find({1, 2, 3})); + Simplex_tree::Filtration_value f013 = st3.filtration(st3.find({0, 1, 3})); + Simplex_tree::Filtration_value f023 = st3.filtration(st3.find({0, 2, 3})); + Simplex_tree::Filtration_value f0123 = st3.filtration(st3.find({0, 1, 2, 3})); + std::cout << "f0123= " << f0123 << " | f012= " << f012 << " - f123= " << f123 << " - f013= " << f013 << + " - f023= " << f023 << std::endl; + BOOST_CHECK(are_almost_the_same(f0123, std::max(f012, std::max(f123, std::max(f013, f023))))); + +} + +using Vector_of_points = std::vector; + +bool is_point_in_list(Vector_of_points points_list, Point point) { + for (auto& point_in_list : points_list) { + if (point_in_list == point) { + return true; // point found + } + } + return false; // point not found +} + +/* Compute the square value of Euclidean distance between two Points given by a range of coordinates. + * The points are assumed to have the same dimension. */ +template< typename Point > +double custom_square_euclidean_distance(const Point &p1,const Point &p2) { + double dist = 0.; + auto it1 = p1.begin(); + auto it2 = p2.begin(); + for (; it1 != p1.end(); ++it1, ++it2) { + double tmp = *it1 - *it2; + dist += tmp*tmp; + } + return dist; +} + +BOOST_AUTO_TEST_CASE(Rips_complex_from_points) { + // ---------------------------------------------------------------------------- + // Init of a list of points + // ---------------------------------------------------------------------------- + Vector_of_points points; + std::vector coords = { 0.0, 0.0, 0.0, 1.0 }; + points.push_back(Point(coords.begin(), coords.end())); + coords = { 0.0, 0.0, 1.0, 0.0 }; + points.push_back(Point(coords.begin(), coords.end())); + coords = { 0.0, 1.0, 0.0, 0.0 }; + points.push_back(Point(coords.begin(), coords.end())); + coords = { 1.0, 0.0, 0.0, 0.0 }; + points.push_back(Point(coords.begin(), coords.end())); + + // ---------------------------------------------------------------------------- + // Init of a rips complex from the list of points + // ---------------------------------------------------------------------------- + Rips_complex rips_complex_from_points(points, 2.0, custom_square_euclidean_distance); + + std::cout << "========== Rips_complex_from_points ==========" << std::endl; + Simplex_tree st; + const int DIMENSION = 3; + BOOST_CHECK(rips_complex_from_points.create_complex(st, DIMENSION)); + + // Another way to check num_simplices + std::cout << "Iterator on rips complex simplices in the filtration order, with [filtration value]:" << std::endl; + int num_simplices = 0; + for (auto f_simplex : st.filtration_simplex_range()) { + num_simplices++; + std::cout << " ( "; + for (auto vertex : st.simplex_vertex_range(f_simplex)) { + std::cout << vertex << " "; + } + std::cout << ") -> " << "[" << st.filtration(f_simplex) << "] "; + std::cout << std::endl; + } + BOOST_CHECK(num_simplices == 15); + std::cout << "st.num_simplices()=" << st.num_simplices() << std::endl; + BOOST_CHECK(st.num_simplices() == 15); + + std::cout << "st.dimension()=" << st.dimension() << std::endl; + BOOST_CHECK(st.dimension() == DIMENSION); + std::cout << "st.num_vertices()=" << st.num_vertices() << std::endl; + BOOST_CHECK(st.num_vertices() == 4); + + for (auto f_simplex : st.filtration_simplex_range()) { + std::cout << "dimension(" << st.dimension(f_simplex) << ") - f = " << st.filtration(f_simplex) << std::endl; + switch (st.dimension(f_simplex)) { + case 0: + BOOST_CHECK(are_almost_the_same(st.filtration(f_simplex), 0.0)); + break; + case 1: + case 2: + case 3: + BOOST_CHECK(are_almost_the_same(st.filtration(f_simplex), 2.0)); + break; + default: + BOOST_CHECK(false); // Shall not happen + break; + } + } +} diff --git a/src/cmake/modules/GUDHI_user_version_target.txt b/src/cmake/modules/GUDHI_user_version_target.txt index 805f0a83..78d3092f 100644 --- a/src/cmake/modules/GUDHI_user_version_target.txt +++ b/src/cmake/modules/GUDHI_user_version_target.txt @@ -48,7 +48,7 @@ if (NOT CMAKE_VERSION VERSION_LESS 2.8.11) add_custom_command(TARGET user_version PRE_BUILD COMMAND ${CMAKE_COMMAND} -E copy_directory ${CMAKE_SOURCE_DIR}/src/GudhUI ${GUDHI_USER_VERSION_DIR}/GudhUI) - set(GUDHI_MODULES "common;Alpha_complex;Bitmap_cubical_complex;Contraction;Hasse_complex;Persistent_cohomology;Simplex_tree;Skeleton_blocker;Witness_complex") + set(GUDHI_MODULES "common;Alpha_complex;Bitmap_cubical_complex;Contraction;Hasse_complex;Persistent_cohomology;Rips_complex;Simplex_tree;Skeleton_blocker;Witness_complex") foreach(GUDHI_MODULE ${GUDHI_MODULES}) # doc files diff --git a/src/common/doc/main_page.h b/src/common/doc/main_page.h index 21cf6925..a38dbcdf 100644 --- a/src/common/doc/main_page.h +++ b/src/common/doc/main_page.h @@ -55,6 +55,23 @@ User manual: \ref cubical_complex - Reference manual: Gudhi::cubical_complex::Bitmap_cubical_complex + \subsection RipsComplexDataStructure Rips complex + \image html "rips_complex_representation.png" "Rips complex representation" + + + + + +
+ Author: Clément Maria
+ Introduced in: GUDHI 1.4.0
+ Copyright: GPL v3
+ 		+ Rips_complex is a simplicial complex constructed from a one skeleton graph.
+ The filtration value of each edge is computed from a user-given distance function and is inserted until a + user-given threshold value.
+ User manual: \ref rips_complex - Reference manual: Gudhi::rips_complex::Rips_complex +
\subsection SimplexTreeDataStructure Simplex tree \image html "Simplex_tree_representation.png" "Simplex tree representation" @@ -288,6 +305,10 @@ make \endverbatim * Persistent_cohomology/periodic_alpha_complex_3d_persistence.cpp * \li * Persistent_cohomology/custom_persistence_sort.cpp + * \li + * Rips_complex/example_one_skeleton_rips_from_points.cpp + * \li + * Rips_complex/example_rips_complex_from_off_file.cpp * * \section Contributions Bug reports and contributions * Please help us improving the quality of the GUDHI library. You may report bugs or suggestions to: @@ -330,6 +351,8 @@ make \endverbatim * @example Persistent_cohomology/rips_multifield_persistence.cpp * @example Persistent_cohomology/rips_persistence.cpp * @example Persistent_cohomology/custom_persistence_sort.cpp + * @example Rips_complex/example_one_skeleton_rips_from_points.cpp + * @example Rips_complex/example_rips_complex_from_off_file.cpp * @example Simplex_tree/mini_simplex_tree.cpp * @example Simplex_tree/simple_simplex_tree.cpp * @example Simplex_tree/simplex_tree_from_alpha_shapes_3.cpp -- cgit v1.2.3 From 29499b02d1b6eafcc6419a0b6b4469152ea20a09 Mon Sep 17 00:00:00 2001 From: vrouvrea Date: Fri, 30 Sep 2016 14:24:41 +0000 Subject: Fix compilation issues git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/rips_complex_module@1596 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: 343308ef914e6a6617890f49a55eb0bb8db15ee9 --- CMakeLists.txt | 1 + src/GudhUI/utils/Persistence_compute.h | 35 ++-- src/Persistent_cohomology/benchmark/CMakeLists.txt | 14 ++ .../benchmark/performance_rips_persistence.cpp | 217 +++++++++++++++++++++ src/Persistent_cohomology/example/CMakeLists.txt | 4 - .../example/alpha_complex_3d_persistence.cpp | 3 +- .../example/alpha_complex_persistence.cpp | 2 + .../example/performance_rips_persistence.cpp | 214 -------------------- .../periodic_alpha_complex_3d_persistence.cpp | 1 + .../persistence_from_simple_simplex_tree.cpp | 17 +- .../example/rips_multifield_persistence.cpp | 21 +- .../example/rips_persistence.cpp | 20 +- .../rips_persistence_via_boundary_matrix.cpp | 50 +++-- ...persistent_cohomology_unit_test_multi_field.cpp | 4 +- src/Rips_complex/include/gudhi/Rips_complex.h | 2 +- src/Simplex_tree/example/simple_simplex_tree.cpp | 13 +- .../example/simplex_tree_from_cliques_of_graph.cpp | 9 +- src/Simplex_tree/include/gudhi/Simplex_tree.h | 2 +- src/Simplex_tree/test/simplex_tree_unit_test.cpp | 69 +++---- .../example/witness_complex_from_file.cpp | 2 +- .../include/gudhi/Witness_complex.h | 2 +- .../test/simple_witness_complex.cpp | 2 +- .../test/witness_complex_points.cpp | 2 +- src/common/include/gudhi/reader_utils.h | 3 +- 24 files changed, 376 insertions(+), 333 deletions(-) create mode 100644 src/Persistent_cohomology/benchmark/CMakeLists.txt create mode 100644 src/Persistent_cohomology/benchmark/performance_rips_persistence.cpp delete mode 100644 src/Persistent_cohomology/example/performance_rips_persistence.cpp (limited to 'CMakeLists.txt') diff --git a/CMakeLists.txt b/CMakeLists.txt index 972df204..cecad3e3 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -134,6 +134,7 @@ else() add_subdirectory(src/Simplex_tree/example) add_subdirectory(src/Persistent_cohomology/test) add_subdirectory(src/Persistent_cohomology/example) + add_subdirectory(src/Persistent_cohomology/benchmark) add_subdirectory(src/Skeleton_blocker/test) add_subdirectory(src/Skeleton_blocker/example) add_subdirectory(src/Contraction/example) diff --git a/src/GudhUI/utils/Persistence_compute.h b/src/GudhUI/utils/Persistence_compute.h index 97165490..fa5bafc1 100644 --- a/src/GudhUI/utils/Persistence_compute.h +++ b/src/GudhUI/utils/Persistence_compute.h @@ -29,6 +29,7 @@ #include #include #include +#include #include @@ -69,22 +70,26 @@ template class Persistence_compute { points.emplace_back(std::move(pt_to_add)); } + using Simplex_tree = Gudhi::Simplex_tree<>; + using Filtration_value = Simplex_tree::Filtration_value; + using Rips_complex = Gudhi::rips_complex::Rips_complex; + using Field_Zp = Gudhi::persistent_cohomology::Field_Zp; + using Persistent_cohomology = Gudhi::persistent_cohomology::Persistent_cohomology; + + Rips_complex rips_complex(points, params.threshold, euclidean_distance); - Graph_t prox_graph = compute_proximity_graph(points, params.threshold, euclidean_distance); - Gudhi::Simplex_tree<> st; - st.insert_graph(prox_graph); - st.expansion(params.max_dim); - - Gudhi::persistent_cohomology::Persistent_cohomology< Gudhi::Simplex_tree<>, - Gudhi::persistent_cohomology::Field_Zp > pcoh(st); - // initializes the coefficient field for homology - pcoh.init_coefficients(params.p); - // put params.min_pers - pcoh.compute_persistent_cohomology(params.min_pers); - stream << "persistence: \n"; - stream << "p dimension birth death: \n"; - - pcoh.output_diagram(stream); + Simplex_tree st; + if (rips_complex.create_complex(st, params.max_dim)) { + Persistent_cohomology pcoh(st); + // initializes the coefficient field for homology + pcoh.init_coefficients(params.p); + // put params.min_pers + pcoh.compute_persistent_cohomology(params.min_pers); + stream << "persistence: \n"; + stream << "p dimension birth death: \n"; + + pcoh.output_diagram(stream); + } } }; diff --git a/src/Persistent_cohomology/benchmark/CMakeLists.txt b/src/Persistent_cohomology/benchmark/CMakeLists.txt new file mode 100644 index 00000000..ea792c89 --- /dev/null +++ b/src/Persistent_cohomology/benchmark/CMakeLists.txt @@ -0,0 +1,14 @@ +cmake_minimum_required(VERSION 2.6) +project(Persistent_cohomology_benchmark) + + +if(GMP_FOUND) + if(GMPXX_FOUND) + add_executable ( performance_rips_persistence EXCLUDE_FROM_ALL performance_rips_persistence.cpp ) + target_link_libraries(performance_rips_persistence ${Boost_SYSTEM_LIBRARY} ${Boost_PROGRAM_OPTIONS_LIBRARY} ${GMPXX_LIBRARIES} ${GMP_LIBRARIES}) + if (TBB_FOUND) + target_link_libraries(performance_rips_persistence ${TBB_LIBRARIES}) + endif(TBB_FOUND) + file(COPY "${CMAKE_SOURCE_DIR}/data/points/Kl.off" DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/) + endif(GMPXX_FOUND) +endif(GMP_FOUND) diff --git a/src/Persistent_cohomology/benchmark/performance_rips_persistence.cpp b/src/Persistent_cohomology/benchmark/performance_rips_persistence.cpp new file mode 100644 index 00000000..99fcad41 --- /dev/null +++ b/src/Persistent_cohomology/benchmark/performance_rips_persistence.cpp @@ -0,0 +1,217 @@ +/* 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 + * + * Copyright (C) 2014 INRIA Sophia Antipolis-Méditerranée (France) + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ + +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +// Types definition +using Simplex_tree = Gudhi::Simplex_tree; +using Filtration_value = Simplex_tree::Filtration_value; +using Rips_complex = Gudhi::rips_complex::Rips_complex; +using Field_Zp = Gudhi::persistent_cohomology::Field_Zp; +using Multi_field = Gudhi::persistent_cohomology::Multi_field; +using Point = std::vector; +using Points_off_reader = Gudhi::Points_off_reader; + +/* Compute the persistent homology of the complex cpx with coefficients in Z/pZ. */ +template< typename FilteredComplex> +void timing_persistence(FilteredComplex & cpx + , int p); + +/* Compute multi-field persistent homology of the complex cpx with coefficients in + * Z/rZ for all prime number r in [p;q].*/ +template< typename FilteredComplex> +void timing_persistence(FilteredComplex & cpx + , int p + , int q); + +/* Timings for the computation of persistent homology with different + * representations of a Rips complex and different coefficient fields. The + * Rips complex is built on a set of 10000 points sampling a Klein bottle embedded + * in dimension 5. + * We represent complexes with a simplex tree and + * with a Hasse diagram. The Hasse diagram represents explicitly all + * codimension 1 incidence relations in the complex, and hence leads to + * a faster computation of persistence because boundaries are precomputed. + * Hovewer, the simplex tree may be constructed directly from a point cloud and + * is more compact. + * We compute persistent homology with coefficient fields Z/2Z and Z/1223Z. + * We present also timings for the computation of multi-field persistent + * homology in all fields Z/rZ for r prime between 2 and 1223. + */ +int main(int argc, char * argv[]) { + std::chrono::time_point start, end; + int elapsed_sec; + { + + std::string off_file_points = "Kl.off"; + Filtration_value threshold = 0.27; + int dim_max = 3; + int p = 2; + int q = 1223; + + // Extract the points from the file off_file_points + Points_off_reader off_reader(off_file_points); + + // Compute the proximity graph of the points + start = std::chrono::system_clock::now(); + Rips_complex rips_complex_from_file(off_reader.get_point_cloud(), threshold, + euclidean_distance); + end = std::chrono::system_clock::now(); + elapsed_sec = std::chrono::duration_cast(end - start).count(); + std::cout << "Compute Rips graph in " << elapsed_sec << " ms.\n"; + + // Construct the Rips complex in a Simplex Tree + Simplex_tree st; + start = std::chrono::system_clock::now(); + + // insert the proximity graph in the simplex tree + // expand the graph until dimension dim_max + rips_complex_from_file.create_complex(st, dim_max); + + end = std::chrono::system_clock::now(); + elapsed_sec = std::chrono::duration_cast(end - start).count(); + std::cout << "Compute Rips complex in " << elapsed_sec << " ms.\n"; + std::cout << " - dimension = " << st.dimension() << std::endl; + std::cout << " - number of simplices = " << st.num_simplices() << std::endl; + + // Sort the simplices in the order of the filtration + start = std::chrono::system_clock::now(); + st.initialize_filtration(); + end = std::chrono::system_clock::now(); + elapsed_sec = std::chrono::duration_cast(end - start).count(); + std::cout << "Order the simplices of the filtration in " << elapsed_sec << " ms.\n"; + + // Copy the keys inside the simplices + start = std::chrono::system_clock::now(); + { + int count = 0; + for (auto sh : st.filtration_simplex_range()) + st.assign_key(sh, count++); + } + end = std::chrono::system_clock::now(); + elapsed_sec = std::chrono::duration_cast(end - start).count(); + std::cout << "Copied the keys inside the simplices in " << elapsed_sec << " ms.\n"; + + // Convert the simplex tree into a hasse diagram + start = std::chrono::system_clock::now(); + Gudhi::Hasse_complex<> hcpx(st); + end = std::chrono::system_clock::now(); + elapsed_sec = std::chrono::duration_cast(end - start).count(); + std::cout << "Convert the simplex tree into a Hasse diagram in " << elapsed_sec << " ms.\n"; + + + std::cout << "Timings when using a simplex tree: \n"; + timing_persistence(st, p); + timing_persistence(st, q); + timing_persistence(st, p, q); + + std::cout << "Timings when using a Hasse complex: \n"; + timing_persistence(hcpx, p); + timing_persistence(hcpx, q); + timing_persistence(hcpx, p, q); + + start = std::chrono::system_clock::now(); + } + end = std::chrono::system_clock::now(); + elapsed_sec = std::chrono::duration_cast(end - start).count(); + std::cout << "Running the complex destructors in " << elapsed_sec << " ms.\n"; + return 0; +} + +template< typename FilteredComplex> +void +timing_persistence(FilteredComplex & cpx + , int p) { + std::chrono::time_point start, end; + int elapsed_sec; + { + start = std::chrono::system_clock::now(); + Gudhi::persistent_cohomology::Persistent_cohomology< FilteredComplex, Field_Zp > pcoh(cpx); + end = std::chrono::system_clock::now(); + elapsed_sec = std::chrono::duration_cast(end - start).count(); + std::cout << " Initialize pcoh in " << elapsed_sec << " ms.\n"; + // initializes the coefficient field for homology + start = std::chrono::system_clock::now(); + pcoh.init_coefficients(p); + end = std::chrono::system_clock::now(); + elapsed_sec = std::chrono::duration_cast(end - start).count(); + std::cout << " Initialize the coefficient field in " << elapsed_sec << " ms.\n"; + + start = std::chrono::system_clock::now(); + + pcoh.compute_persistent_cohomology(INFINITY); + + end = std::chrono::system_clock::now(); + elapsed_sec = std::chrono::duration_cast(end - start).count(); + std::cout << " Compute persistent homology in Z/" << p << "Z in " << elapsed_sec << " ms.\n"; + start = std::chrono::system_clock::now(); + } + end = std::chrono::system_clock::now(); + elapsed_sec = std::chrono::duration_cast(end - start).count(); + std::cout << " Run the persistence destructors in " << elapsed_sec << " ms.\n"; +} + +template< typename FilteredComplex> +void +timing_persistence(FilteredComplex & cpx + , int p + , int q) { + std::chrono::time_point start, end; + int elapsed_sec; + { + start = std::chrono::system_clock::now(); + Gudhi::persistent_cohomology::Persistent_cohomology< FilteredComplex, Multi_field > pcoh(cpx); + end = std::chrono::system_clock::now(); + elapsed_sec = std::chrono::duration_cast(end - start).count(); + std::cout << " Initialize pcoh in " << elapsed_sec << " ms.\n"; + // initializes the coefficient field for homology + start = std::chrono::system_clock::now(); + pcoh.init_coefficients(p, q); + end = std::chrono::system_clock::now(); + elapsed_sec = std::chrono::duration_cast(end - start).count(); + std::cout << " Initialize the coefficient field in " << elapsed_sec << " ms.\n"; + // compute persistent homology, disgarding persistent features of life shorter than min_persistence + + start = std::chrono::system_clock::now(); + + pcoh.compute_persistent_cohomology(INFINITY); + + end = std::chrono::system_clock::now(); + elapsed_sec = std::chrono::duration_cast(end - start).count(); + std::cout << " Compute multi-field persistent homology in all coefficient fields Z/pZ " + << "with p in [" << p << ";" << q << "] in " << elapsed_sec << " ms.\n"; + start = std::chrono::system_clock::now(); + } + end = std::chrono::system_clock::now(); + elapsed_sec = std::chrono::duration_cast(end - start).count(); + std::cout << " Run the persistence destructors in " << elapsed_sec << " ms.\n"; +} diff --git a/src/Persistent_cohomology/example/CMakeLists.txt b/src/Persistent_cohomology/example/CMakeLists.txt index d97d1b63..75f3c5d4 100644 --- a/src/Persistent_cohomology/example/CMakeLists.txt +++ b/src/Persistent_cohomology/example/CMakeLists.txt @@ -39,13 +39,9 @@ if(GMP_FOUND) if(GMPXX_FOUND) add_executable(rips_multifield_persistence rips_multifield_persistence.cpp ) target_link_libraries(rips_multifield_persistence ${Boost_SYSTEM_LIBRARY} ${Boost_PROGRAM_OPTIONS_LIBRARY} ${GMPXX_LIBRARIES} ${GMP_LIBRARIES}) - add_executable ( performance_rips_persistence performance_rips_persistence.cpp ) - target_link_libraries(performance_rips_persistence ${Boost_SYSTEM_LIBRARY} ${Boost_PROGRAM_OPTIONS_LIBRARY} ${GMPXX_LIBRARIES} ${GMP_LIBRARIES}) if (TBB_FOUND) target_link_libraries(rips_multifield_persistence ${TBB_LIBRARIES}) - target_link_libraries(performance_rips_persistence ${TBB_LIBRARIES}) endif(TBB_FOUND) - add_test(rips_multifield_persistence_2_71 ${CMAKE_CURRENT_BINARY_DIR}/rips_multifield_persistence ${CMAKE_SOURCE_DIR}/data/points/Kl.txt -r 0.2 -d 3 -p 2 -q 71 -m 100) endif(GMPXX_FOUND) endif(GMP_FOUND) diff --git a/src/Persistent_cohomology/example/alpha_complex_3d_persistence.cpp b/src/Persistent_cohomology/example/alpha_complex_3d_persistence.cpp index 48fbb91a..f7e8f800 100644 --- a/src/Persistent_cohomology/example/alpha_complex_3d_persistence.cpp +++ b/src/Persistent_cohomology/example/alpha_complex_3d_persistence.cpp @@ -64,6 +64,7 @@ typedef std::list Vertex_list; // gudhi type definition typedef Gudhi::Simplex_tree ST; +typedef ST::Filtration_value Filtration_value; typedef ST::Vertex_handle Simplex_tree_vertex; typedef std::map Alpha_shape_simplex_tree_map; typedef std::pair Alpha_shape_simplex_tree_pair; @@ -132,7 +133,7 @@ int main(int argc, char * const argv[]) { int coeff_field_characteristic = atoi(argv[2]); Filtration_value min_persistence = 0.0; - int returnedScanValue = sscanf(argv[3], "%lf", &min_persistence); + int returnedScanValue = sscanf(argv[3], "%f", &min_persistence); if ((returnedScanValue == EOF) || (min_persistence < -1.0)) { std::cerr << "Error: " << argv[3] << " is not correct\n"; usage(argv[0]); diff --git a/src/Persistent_cohomology/example/alpha_complex_persistence.cpp b/src/Persistent_cohomology/example/alpha_complex_persistence.cpp index cb181936..24315144 100644 --- a/src/Persistent_cohomology/example/alpha_complex_persistence.cpp +++ b/src/Persistent_cohomology/example/alpha_complex_persistence.cpp @@ -9,6 +9,8 @@ #include #include // for numeric_limits +using Filtration_value = double; + void program_options(int argc, char * argv[] , std::string & off_file_points , std::string & output_file_diag diff --git a/src/Persistent_cohomology/example/performance_rips_persistence.cpp b/src/Persistent_cohomology/example/performance_rips_persistence.cpp deleted file mode 100644 index b4d282ac..00000000 --- a/src/Persistent_cohomology/example/performance_rips_persistence.cpp +++ /dev/null @@ -1,214 +0,0 @@ -/* 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 - * - * Copyright (C) 2014 INRIA Sophia Antipolis-Méditerranée (France) - * - * This program is free software: you can redistribute it and/or modify - * it under the terms of the GNU General Public License as published by - * the Free Software Foundation, either version 3 of the License, or - * (at your option) any later version. - * - * This program is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - * GNU General Public License for more details. - * - * You should have received a copy of the GNU General Public License - * along with this program. If not, see . - */ - -#include -#include -#include -#include -#include -#include -#include - -#include -#include -#include - -using namespace Gudhi; -using namespace Gudhi::persistent_cohomology; - -/* Compute the persistent homology of the complex cpx with coefficients in Z/pZ. */ -template< typename FilteredComplex> -void timing_persistence(FilteredComplex & cpx - , int p); - -/* Compute multi-field persistent homology of the complex cpx with coefficients in - * Z/rZ for all prime number r in [p;q].*/ -template< typename FilteredComplex> -void timing_persistence(FilteredComplex & cpx - , int p - , int q); - -/* Timings for the computation of persistent homology with different - * representations of a Rips complex and different coefficient fields. The - * Rips complex is built on a set of 10000 points sampling a Klein bottle embedded - * in dimension 5. - * We represent complexes with a simplex tree and - * with a Hasse diagram. The Hasse diagram represents explicitly all - * codimension 1 incidence relations in the complex, and hence leads to - * a faster computation of persistence because boundaries are precomputed. - * Hovewer, the simplex tree may be constructed directly from a point cloud and - * is more compact. - * We compute persistent homology with coefficient fields Z/2Z and Z/1223Z. - * We present also timings for the computation of multi-field persistent - * homology in all fields Z/rZ for r prime between 2 and 1223. - */ -int main(int argc, char * argv[]) { - std::chrono::time_point start, end; - int elapsed_sec; - { - - std::string filepoints = "../../../data/points/Kl.txt"; - Filtration_value threshold = 0.27; - int dim_max = 3; - int p = 2; - int q = 1223; - - // Extract the points from the file filepoints - typedef std::vector Point_t; - std::vector< Point_t > points; - read_points(filepoints, points); - - // Compute the proximity graph of the points - start = std::chrono::system_clock::now(); - Graph_t prox_graph = compute_proximity_graph(points, threshold - , euclidean_distance); - end = std::chrono::system_clock::now(); - elapsed_sec = std::chrono::duration_cast(end - start).count(); - std::cout << "Compute Rips graph in " << elapsed_sec << " ms.\n"; - - // Construct the Rips complex in a Simplex Tree - Simplex_tree st; - start = std::chrono::system_clock::now(); - - // insert the proximity graph in the simplex tree - st.insert_graph(prox_graph); - // expand the graph until dimension dim_max - st.expansion(dim_max); - - end = std::chrono::system_clock::now(); - elapsed_sec = std::chrono::duration_cast(end - start).count(); - std::cout << "Compute Rips complex in " << elapsed_sec << " ms.\n"; - std::cout << " - dimension = " << st.dimension() << std::endl; - std::cout << " - number of simplices = " << st.num_simplices() << std::endl; - - // Sort the simplices in the order of the filtration - start = std::chrono::system_clock::now(); - st.initialize_filtration(); - end = std::chrono::system_clock::now(); - elapsed_sec = std::chrono::duration_cast(end - start).count(); - std::cout << "Order the simplices of the filtration in " << elapsed_sec << " ms.\n"; - - // Copy the keys inside the simplices - start = std::chrono::system_clock::now(); - { - int count = 0; - for (auto sh : st.filtration_simplex_range()) - st.assign_key(sh, count++); - } - end = std::chrono::system_clock::now(); - elapsed_sec = std::chrono::duration_cast(end - start).count(); - std::cout << "Copied the keys inside the simplices in " << elapsed_sec << " ms.\n"; - - // Convert the simplex tree into a hasse diagram - start = std::chrono::system_clock::now(); - Hasse_complex<> hcpx(st); - end = std::chrono::system_clock::now(); - elapsed_sec = std::chrono::duration_cast(end - start).count(); - std::cout << "Convert the simplex tree into a Hasse diagram in " << elapsed_sec << " ms.\n"; - - - std::cout << "Timings when using a simplex tree: \n"; - timing_persistence(st, p); - timing_persistence(st, q); - timing_persistence(st, p, q); - - std::cout << "Timings when using a Hasse complex: \n"; - timing_persistence(hcpx, p); - timing_persistence(hcpx, q); - timing_persistence(hcpx, p, q); - - start = std::chrono::system_clock::now(); - } - end = std::chrono::system_clock::now(); - elapsed_sec = std::chrono::duration_cast(end - start).count(); - std::cout << "Running the complex destructors in " << elapsed_sec << " ms.\n"; - return 0; -} - -template< typename FilteredComplex> -void -timing_persistence(FilteredComplex & cpx - , int p) { - std::chrono::time_point start, end; - int elapsed_sec; - { - start = std::chrono::system_clock::now(); - Persistent_cohomology< FilteredComplex, Field_Zp > pcoh(cpx); - end = std::chrono::system_clock::now(); - elapsed_sec = std::chrono::duration_cast(end - start).count(); - std::cout << " Initialize pcoh in " << elapsed_sec << " ms.\n"; - // initializes the coefficient field for homology - start = std::chrono::system_clock::now(); - pcoh.init_coefficients(p); - end = std::chrono::system_clock::now(); - elapsed_sec = std::chrono::duration_cast(end - start).count(); - std::cout << " Initialize the coefficient field in " << elapsed_sec << " ms.\n"; - - start = std::chrono::system_clock::now(); - - pcoh.compute_persistent_cohomology(INFINITY); - - end = std::chrono::system_clock::now(); - elapsed_sec = std::chrono::duration_cast(end - start).count(); - std::cout << " Compute persistent homology in Z/" << p << "Z in " << elapsed_sec << " ms.\n"; - start = std::chrono::system_clock::now(); - } - end = std::chrono::system_clock::now(); - elapsed_sec = std::chrono::duration_cast(end - start).count(); - std::cout << " Run the persistence destructors in " << elapsed_sec << " ms.\n"; -} - -template< typename FilteredComplex> -void -timing_persistence(FilteredComplex & cpx - , int p - , int q) { - std::chrono::time_point start, end; - int elapsed_sec; - { - start = std::chrono::system_clock::now(); - Persistent_cohomology< FilteredComplex, Multi_field > pcoh(cpx); - end = std::chrono::system_clock::now(); - elapsed_sec = std::chrono::duration_cast(end - start).count(); - std::cout << " Initialize pcoh in " << elapsed_sec << " ms.\n"; - // initializes the coefficient field for homology - start = std::chrono::system_clock::now(); - pcoh.init_coefficients(p, q); - end = std::chrono::system_clock::now(); - elapsed_sec = std::chrono::duration_cast(end - start).count(); - std::cout << " Initialize the coefficient field in " << elapsed_sec << " ms.\n"; - // compute persistent homology, disgarding persistent features of life shorter than min_persistence - - start = std::chrono::system_clock::now(); - - pcoh.compute_persistent_cohomology(INFINITY); - - end = std::chrono::system_clock::now(); - elapsed_sec = std::chrono::duration_cast(end - start).count(); - std::cout << " Compute multi-field persistent homology in all coefficient fields Z/pZ " - << "with p in [" << p << ";" << q << "] in " << elapsed_sec << " ms.\n"; - start = std::chrono::system_clock::now(); - } - end = std::chrono::system_clock::now(); - elapsed_sec = std::chrono::duration_cast(end - start).count(); - std::cout << " Run the persistence destructors in " << elapsed_sec << " ms.\n"; -} diff --git a/src/Persistent_cohomology/example/periodic_alpha_complex_3d_persistence.cpp b/src/Persistent_cohomology/example/periodic_alpha_complex_3d_persistence.cpp index a199fea1..5184ef52 100644 --- a/src/Persistent_cohomology/example/periodic_alpha_complex_3d_persistence.cpp +++ b/src/Persistent_cohomology/example/periodic_alpha_complex_3d_persistence.cpp @@ -70,6 +70,7 @@ using Vertex_list = std::list; // gudhi type definition using ST = Gudhi::Simplex_tree; +using Filtration_value = ST::Filtration_value; using Simplex_tree_vertex = ST::Vertex_handle; using Alpha_shape_simplex_tree_map = std::map; using Alpha_shape_simplex_tree_pair = std::pair; diff --git a/src/Persistent_cohomology/example/persistence_from_simple_simplex_tree.cpp b/src/Persistent_cohomology/example/persistence_from_simple_simplex_tree.cpp index ba772f04..817aac4d 100644 --- a/src/Persistent_cohomology/example/persistence_from_simple_simplex_tree.cpp +++ b/src/Persistent_cohomology/example/persistence_from_simple_simplex_tree.cpp @@ -29,13 +29,12 @@ #include #include -using namespace Gudhi; -using namespace Gudhi::persistent_cohomology; - -typedef std::vector< Vertex_handle > typeVectorVertex; -typedef std::pair typeSimplex; -typedef std::pair< Simplex_tree<>::Simplex_handle, bool > typePairSimplexBool; -typedef Simplex_tree<> typeST; +// Types definition +using Simplex_tree = Gudhi::Simplex_tree; +using Filtration_value = Simplex_tree::Filtration_value; +using Field_Zp = Gudhi::persistent_cohomology::Field_Zp; +using Persistent_cohomology = Gudhi::persistent_cohomology::Persistent_cohomology; +using typeVectorVertex = std::vector< Simplex_tree::Vertex_handle >; void usage(char * const progName) { std::cerr << "Usage: " << progName << " coeff_field_characteristic[integer > 0] min_persistence[float >= -1.0]\n"; @@ -66,7 +65,7 @@ int main(int argc, char * const argv[]) { // TEST OF INSERTION std::cout << "********************************************************************" << std::endl; std::cout << "TEST OF INSERTION" << std::endl; - typeST st; + Simplex_tree st; // ++ FIRST std::cout << " - INSERT (0,1,2)" << std::endl; @@ -166,7 +165,7 @@ int main(int argc, char * const argv[]) { std::cout << "**************************************************************" << std::endl; // Compute the persistence diagram of the complex - persistent_cohomology::Persistent_cohomology< Simplex_tree<>, Field_Zp > pcoh(st); + Persistent_cohomology pcoh(st); // initializes the coefficient field for homology pcoh.init_coefficients(coeff_field_characteristic); diff --git a/src/Persistent_cohomology/example/rips_multifield_persistence.cpp b/src/Persistent_cohomology/example/rips_multifield_persistence.cpp index d4a07bb6..f4adc7a9 100644 --- a/src/Persistent_cohomology/example/rips_multifield_persistence.cpp +++ b/src/Persistent_cohomology/example/rips_multifield_persistence.cpp @@ -25,14 +25,21 @@ #include #include #include +#include #include #include #include - -typedef double Filtration_value; +// Types definition +using Simplex_tree = Gudhi::Simplex_tree; +using Filtration_value = Simplex_tree::Filtration_value; +using Rips_complex = Gudhi::rips_complex::Rips_complex; +using Multi_field = Gudhi::persistent_cohomology::Multi_field; +using Persistent_cohomology = Gudhi::persistent_cohomology::Persistent_cohomology; +using Point = std::vector; +using Points_off_reader = Gudhi::Points_off_reader; void program_options(int argc, char * argv[] , std::string & off_file_points @@ -54,22 +61,22 @@ int main(int argc, char * argv[]) { program_options(argc, argv, off_file_points, filediag, threshold, dim_max, min_p, max_p, min_persistence); - Gudhi::rips_complex::Rips_complex<> rips_complex_from_file(off_file_points); + Points_off_reader off_reader(off_file_points); + Rips_complex rips_complex_from_file(off_reader.get_point_cloud(), threshold, + euclidean_distance); // Construct the Rips complex in a Simplex Tree - using Simplex_tree = Gudhi::Simplex_tree; Simplex_tree simplex_tree; - if (rips_complex_from_file.create_complex(simplex_tree, threshold, dim_max, euclidean_distance>)) { + if (rips_complex_from_file.create_complex(simplex_tree, dim_max)) { std::cout << "The complex contains " << simplex_tree.num_simplices() << " simplices \n"; std::cout << " and has dimension " << simplex_tree.dimension() << " \n"; // Sort the simplices in the order of the filtration simplex_tree.initialize_filtration(); - using Multi_field = Gudhi::persistent_cohomology::Multi_field; // Compute the persistence diagram of the complex - Gudhi::persistent_cohomology::Persistent_cohomology pcoh(simplex_tree); + Persistent_cohomology pcoh(simplex_tree); // initializes the coefficient field for homology pcoh.init_coefficients(min_p, max_p); diff --git a/src/Persistent_cohomology/example/rips_persistence.cpp b/src/Persistent_cohomology/example/rips_persistence.cpp index e3ab5927..97bab14c 100644 --- a/src/Persistent_cohomology/example/rips_persistence.cpp +++ b/src/Persistent_cohomology/example/rips_persistence.cpp @@ -25,6 +25,7 @@ #include #include #include +#include #include @@ -32,7 +33,14 @@ #include #include // infinity -typedef double Filtration_value; +// Types definition +using Simplex_tree = Gudhi::Simplex_tree; +using Filtration_value = Simplex_tree::Filtration_value; +using Rips_complex = Gudhi::rips_complex::Rips_complex; +using Field_Zp = Gudhi::persistent_cohomology::Field_Zp; +using Persistent_cohomology = Gudhi::persistent_cohomology::Persistent_cohomology; +using Point = std::vector; +using Points_off_reader = Gudhi::Points_off_reader; void program_options(int argc, char * argv[] , std::string & off_file_points @@ -52,22 +60,22 @@ int main(int argc, char * argv[]) { program_options(argc, argv, off_file_points, filediag, threshold, dim_max, p, min_persistence); - Gudhi::rips_complex::Rips_complex<> rips_complex_from_file(off_file_points); + Points_off_reader off_reader(off_file_points); + Rips_complex rips_complex_from_file(off_reader.get_point_cloud(), threshold, + euclidean_distance); // Construct the Rips complex in a Simplex Tree - using Simplex_tree = Gudhi::Simplex_tree; Simplex_tree simplex_tree; - if (rips_complex_from_file.create_complex(simplex_tree, threshold, dim_max, euclidean_distance>)) { + if (rips_complex_from_file.create_complex(simplex_tree, dim_max)) { std::cout << "The complex contains " << simplex_tree.num_simplices() << " simplices \n"; std::cout << " and has dimension " << simplex_tree.dimension() << " \n"; // Sort the simplices in the order of the filtration simplex_tree.initialize_filtration(); - using Field_Zp = Gudhi::persistent_cohomology::Field_Zp; // Compute the persistence diagram of the complex - Gudhi::persistent_cohomology::Persistent_cohomology pcoh(simplex_tree); + Persistent_cohomology pcoh(simplex_tree); // initializes the coefficient field for homology pcoh.init_coefficients(p); diff --git a/src/Persistent_cohomology/example/rips_persistence_via_boundary_matrix.cpp b/src/Persistent_cohomology/example/rips_persistence_via_boundary_matrix.cpp index 4c6656f5..991ed4f0 100644 --- a/src/Persistent_cohomology/example/rips_persistence_via_boundary_matrix.cpp +++ b/src/Persistent_cohomology/example/rips_persistence_via_boundary_matrix.cpp @@ -21,12 +21,12 @@ * along with this program. If not, see . */ -#include -#include -#include #include #include +#include #include +#include +#include #include @@ -44,14 +44,16 @@ // // //////////////////////////////////////////////////////////////// -using namespace Gudhi; -using namespace Gudhi::persistent_cohomology; - -typedef int Vertex_handle; -typedef double Filtration_value; +// Types definition +using Simplex_tree = Gudhi::Simplex_tree<>; +using Filtration_value = Simplex_tree::Filtration_value; +using Rips_complex = Gudhi::rips_complex::Rips_complex; +using Field_Zp = Gudhi::persistent_cohomology::Field_Zp; +using Point = std::vector; +using Points_off_reader = Gudhi::Points_off_reader; void program_options(int argc, char * argv[] - , std::string & filepoints + , std::string & off_file_points , std::string & filediag , Filtration_value & threshold , int & dim_max @@ -59,30 +61,22 @@ void program_options(int argc, char * argv[] , Filtration_value & min_persistence); int main(int argc, char * argv[]) { - std::string filepoints; + std::string off_file_points; std::string filediag; Filtration_value threshold; int dim_max; int p; Filtration_value min_persistence; - program_options(argc, argv, filepoints, filediag, threshold, dim_max, p, min_persistence); - - // Extract the points from the file filepoints - typedef std::vector Point_t; - std::vector< Point_t > points; - read_points(filepoints, points); + program_options(argc, argv, off_file_points, filediag, threshold, dim_max, p, min_persistence); - // Compute the proximity graph of the points - Graph_t prox_graph = compute_proximity_graph(points, threshold - , euclidean_distance); + Points_off_reader off_reader(off_file_points); + Rips_complex rips_complex_from_file(off_reader.get_point_cloud(), threshold, + euclidean_distance); // Construct the Rips complex in a Simplex Tree - Simplex_tree<>& st = *new Simplex_tree<>; - // insert the proximity graph in the simplex tree - st.insert_graph(prox_graph); - // expand the graph until dimension dim_max - st.expansion(dim_max); + Simplex_tree& st = *new Simplex_tree; + rips_complex_from_file.create_complex(st, dim_max); std::cout << "The complex contains " << st.num_simplices() << " simplices \n"; std::cout << " and has dimension " << st.dimension() << " \n"; @@ -99,7 +93,7 @@ int main(int argc, char * argv[]) { st.assign_key(sh, count++); // Convert to a more convenient representation. - Hasse_complex<> hcpx(st); + Gudhi::Hasse_complex<> hcpx(st); #ifdef GUDHI_USE_TBB ts.terminate(); @@ -109,7 +103,7 @@ int main(int argc, char * argv[]) { delete &st; // Compute the persistence diagram of the complex - persistent_cohomology::Persistent_cohomology< Hasse_complex<>, Field_Zp > pcoh(hcpx); + Gudhi::persistent_cohomology::Persistent_cohomology< Gudhi::Hasse_complex<>, Field_Zp > pcoh(hcpx); // initializes the coefficient field for homology pcoh.init_coefficients(p); @@ -126,7 +120,7 @@ int main(int argc, char * argv[]) { } void program_options(int argc, char * argv[] - , std::string & filepoints + , std::string & off_file_points , std::string & filediag , Filtration_value & threshold , int & dim_max @@ -135,7 +129,7 @@ void program_options(int argc, char * argv[] namespace po = boost::program_options; po::options_description hidden("Hidden options"); hidden.add_options() - ("input-file", po::value(&filepoints), + ("input-file", po::value(&off_file_points), "Name of file containing a point set. Format is one point per line: X1 ... Xd "); po::options_description visible("Allowed options", 100); diff --git a/src/Persistent_cohomology/test/persistent_cohomology_unit_test_multi_field.cpp b/src/Persistent_cohomology/test/persistent_cohomology_unit_test_multi_field.cpp index 703682e1..1a6e3296 100644 --- a/src/Persistent_cohomology/test/persistent_cohomology_unit_test_multi_field.cpp +++ b/src/Persistent_cohomology/test/persistent_cohomology_unit_test_multi_field.cpp @@ -21,7 +21,7 @@ using namespace boost::unit_test; typedef Simplex_tree<> typeST; -std::string test_rips_persistence(int min_coefficient, int max_coefficient, int min_persistence) { +std::string test_rips_persistence(int min_coefficient, int max_coefficient, double min_persistence) { // file name is given as parameter from CMakeLists.txt const std::string inputFile(framework::master_test_suite().argv[1]); @@ -74,7 +74,7 @@ void test_rips_persistence_in_dimension(int min_dimension, int max_dimension) { std::cout << "********************************************************************" << std::endl; std::cout << "TEST OF RIPS_PERSISTENT_COHOMOLOGY_MULTI_FIELD MIN_DIM=" << min_dimension << " MAX_DIM=" << max_dimension << " MIN_PERS=0" << std::endl; - std::string str_rips_persistence = test_rips_persistence(min_dimension, max_dimension, static_cast (0.0)); + std::string str_rips_persistence = test_rips_persistence(min_dimension, max_dimension, 0.0); std::cout << "str_rips_persistence=" << str_rips_persistence << std::endl; BOOST_CHECK(str_rips_persistence.find(value0) != std::string::npos); // Check found diff --git a/src/Rips_complex/include/gudhi/Rips_complex.h b/src/Rips_complex/include/gudhi/Rips_complex.h index 10f674e5..9eb22d03 100644 --- a/src/Rips_complex/include/gudhi/Rips_complex.h +++ b/src/Rips_complex/include/gudhi/Rips_complex.h @@ -47,7 +47,7 @@ namespace rips_complex { * \ingroup rips_complex * * \details - * The data structure is a 1-skeleton graph constructed from a point cloud, containing edges when the edge length is + * The data structure is a one skeleton graph constructed from a point cloud, containing edges when the edge length is * less or equal to a given threshold. Edge length is computed from a user given function. * * The complex is a template class requiring a Filtration_value type. diff --git a/src/Simplex_tree/example/simple_simplex_tree.cpp b/src/Simplex_tree/example/simple_simplex_tree.cpp index 5146b906..bf6dc470 100644 --- a/src/Simplex_tree/example/simple_simplex_tree.cpp +++ b/src/Simplex_tree/example/simple_simplex_tree.cpp @@ -27,10 +27,11 @@ #include // for pair #include -using namespace Gudhi; - -typedef std::vector< Vertex_handle > typeVectorVertex; -typedef std::pair< Simplex_tree<>::Simplex_handle, bool > typePairSimplexBool; +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 >; int main(int argc, char * const argv[]) { const Filtration_value FIRST_FILTRATION_VALUE = 0.1; @@ -42,7 +43,7 @@ int main(int argc, char * const argv[]) { std::cout << "********************************************************************" << std::endl; std::cout << "EXAMPLE OF SIMPLE INSERTION" << std::endl; // Construct the Simplex Tree - Simplex_tree<> simplexTree; + Simplex_tree simplexTree; /* Simplex to be inserted: */ /* 1 */ @@ -212,7 +213,7 @@ int main(int argc, char * const argv[]) { // ------------------------------------------------------------------------------------------------------------------ // Find in the simplex_tree // ------------------------------------------------------------------------------------------------------------------ - Simplex_tree<>::Simplex_handle simplexFound = simplexTree.find(secondSimplexVector); + Simplex_tree::Simplex_handle simplexFound = simplexTree.find(secondSimplexVector); std::cout << "**************IS THE SIMPLEX {1} IN THE SIMPLEX TREE ?\n"; if (simplexFound != simplexTree.null_simplex()) std::cout << "***+ YES IT IS!\n"; diff --git a/src/Simplex_tree/example/simplex_tree_from_cliques_of_graph.cpp b/src/Simplex_tree/example/simplex_tree_from_cliques_of_graph.cpp index 58085014..8d729c56 100644 --- a/src/Simplex_tree/example/simplex_tree_from_cliques_of_graph.cpp +++ b/src/Simplex_tree/example/simplex_tree_from_cliques_of_graph.cpp @@ -29,6 +29,13 @@ using namespace Gudhi; +typedef int Vertex_handle; +typedef double Filtration_value; +typedef boost::adjacency_list < boost::vecS, boost::vecS, boost::undirectedS, + boost::property < vertex_filtration_t, Filtration_value >, + boost::property < edge_filtration_t, Filtration_value > > Graph_t; +typedef std::pair< Vertex_handle, Vertex_handle > Edge_t; + int main(int argc, char * const argv[]) { if (argc != 3) { std::cerr << "Usage: " << argv[0] @@ -43,7 +50,7 @@ int main(int argc, char * const argv[]) { Simplex_tree<> st; start = clock(); - auto g = read_graph(filegraph); + auto g = read_graph(filegraph); // insert the graph in the simplex tree as 1-skeleton st.insert_graph(g); end = clock(); diff --git a/src/Simplex_tree/include/gudhi/Simplex_tree.h b/src/Simplex_tree/include/gudhi/Simplex_tree.h index 206b2fba..63e3f0e5 100644 --- a/src/Simplex_tree/include/gudhi/Simplex_tree.h +++ b/src/Simplex_tree/include/gudhi/Simplex_tree.h @@ -28,7 +28,7 @@ #include #include -//#include +#include #include #include diff --git a/src/Simplex_tree/test/simplex_tree_unit_test.cpp b/src/Simplex_tree/test/simplex_tree_unit_test.cpp index 28bf202b..b06d7ec9 100644 --- a/src/Simplex_tree/test/simplex_tree_unit_test.cpp +++ b/src/Simplex_tree/test/simplex_tree_unit_test.cpp @@ -1,4 +1,5 @@ #include +#include #include #include #include // std::pair, std::make_pair @@ -19,19 +20,19 @@ using namespace Gudhi; typedef boost::mpl::list, Simplex_tree> list_of_tested_variants; -const Vertex_handle DEFAULT_VERTEX_HANDLE = (const Vertex_handle) - 1; -const Filtration_value DEFAULT_FILTRATION_VALUE = (const Filtration_value) 0.0; template void test_empty_simplex_tree(typeST& tst) { - BOOST_CHECK(tst.null_vertex() == DEFAULT_VERTEX_HANDLE); - BOOST_CHECK(tst.filtration() == DEFAULT_FILTRATION_VALUE); + typedef typename typeST::Vertex_handle Vertex_handle; + const Vertex_handle DEFAULT_VERTEX_VALUE = Vertex_handle(- 1); + BOOST_CHECK(tst.null_vertex() == DEFAULT_VERTEX_VALUE); + BOOST_CHECK(tst.filtration() == 0.0); BOOST_CHECK(tst.num_vertices() == (size_t) 0); BOOST_CHECK(tst.num_simplices() == (size_t) 0); typename typeST::Siblings* STRoot = tst.root(); BOOST_CHECK(STRoot != nullptr); BOOST_CHECK(STRoot->oncles() == nullptr); - BOOST_CHECK(STRoot->parent() == DEFAULT_VERTEX_HANDLE); + BOOST_CHECK(STRoot->parent() == DEFAULT_VERTEX_VALUE); BOOST_CHECK(tst.dimension() == -1); } @@ -59,7 +60,7 @@ void test_iterators_on_empty_simplex_tree(typeST& tst) { BOOST_AUTO_TEST_CASE_TEMPLATE(simplex_tree_when_empty, typeST, list_of_tested_variants) { typedef std::pair typePairSimplexBool; - typedef std::vector typeVectorVertex; + typedef std::vector typeVectorVertex; std::cout << "********************************************************************" << std::endl; std::cout << "TEST OF DEFAULT CONSTRUCTOR" << std::endl; @@ -72,8 +73,7 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(simplex_tree_when_empty, typeST, list_of_tested_va std::cout << "TEST OF EMPTY INSERTION" << std::endl; typeVectorVertex simplexVectorEmpty; BOOST_CHECK(simplexVectorEmpty.empty() == true); - typePairSimplexBool returnEmptyValue = st.insert_simplex(simplexVectorEmpty, - DEFAULT_FILTRATION_VALUE); + typePairSimplexBool returnEmptyValue = st.insert_simplex(simplexVectorEmpty, 0.0); BOOST_CHECK(returnEmptyValue.first == typename typeST::Simplex_handle(nullptr)); BOOST_CHECK(returnEmptyValue.second == true); @@ -141,12 +141,13 @@ void test_simplex_tree_contains(typeST& simplexTree, typeSimplex& simplex, int p template void test_simplex_tree_insert_returns_true(const typePairSimplexBool& returnValue) { BOOST_CHECK(returnValue.second == true); - typename typeST::Simplex_handle shReturned = returnValue.first; // Simplex_handle = boost::container::flat_map< Vertex_handle, Node >::iterator + // Simplex_handle = boost::container::flat_map< typeST::Vertex_handle, Node >::iterator + typename typeST::Simplex_handle shReturned = returnValue.first; BOOST_CHECK(shReturned != typename typeST::Simplex_handle(nullptr)); } // Global variables -Filtration_value max_fil = DEFAULT_FILTRATION_VALUE; +double max_fil = 0.0; int dim_max = -1; template @@ -179,8 +180,9 @@ void set_and_test_simplex_tree_dim_fil(typeST& simplexTree, int vectorSize, cons } BOOST_AUTO_TEST_CASE_TEMPLATE(simplex_tree_insertion, typeST, list_of_tested_variants) { + typedef typename typeST::Filtration_value Filtration_value; typedef std::pair typePairSimplexBool; - typedef std::vector typeVectorVertex; + typedef std::vector typeVectorVertex; typedef std::pair typeSimplex; const Filtration_value FIRST_FILTRATION_VALUE = 0.1; const Filtration_value SECOND_FILTRATION_VALUE = 0.2; @@ -188,7 +190,7 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(simplex_tree_insertion, typeST, list_of_tested_var const Filtration_value FOURTH_FILTRATION_VALUE = 0.4; // reset since we run the test several times dim_max = -1; - max_fil = DEFAULT_FILTRATION_VALUE; + max_fil = 0.0; // TEST OF INSERTION std::cout << "********************************************************************" << std::endl; @@ -303,7 +305,8 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(simplex_tree_insertion, typeST, list_of_tested_var returnValue = st.insert_simplex(tenthSimplex.first, tenthSimplex.second); BOOST_CHECK(returnValue.second == false); - typename typeST::Simplex_handle shReturned = returnValue.first; // Simplex_handle = boost::container::flat_map< Vertex_handle, Node >::iterator + // Simplex_handle = boost::container::flat_map< typeST::Vertex_handle, Node >::iterator + typename typeST::Simplex_handle shReturned = returnValue.first; BOOST_CHECK(shReturned == typename typeST::Simplex_handle(nullptr)); BOOST_CHECK(st.num_vertices() == (size_t) 4); // Not incremented !! BOOST_CHECK(st.dimension() == dim_max); @@ -317,7 +320,8 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(simplex_tree_insertion, typeST, list_of_tested_var returnValue = st.insert_simplex(eleventhSimplex.first, eleventhSimplex.second); BOOST_CHECK(returnValue.second == false); - shReturned = returnValue.first; // Simplex_handle = boost::container::flat_map< Vertex_handle, Node >::iterator + // Simplex_handle = boost::container::flat_map< typeST::Vertex_handle, Node >::iterator + shReturned = returnValue.first; BOOST_CHECK(shReturned == typename typeST::Simplex_handle(nullptr)); BOOST_CHECK(st.num_vertices() == (size_t) 4); // Not incremented !! BOOST_CHECK(st.dimension() == dim_max); @@ -375,8 +379,8 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(simplex_tree_insertion, typeST, list_of_tested_var BOOST_AUTO_TEST_CASE_TEMPLATE(NSimplexAndSubfaces_tree_insertion, typeST, list_of_tested_variants) { typedef std::pair typePairSimplexBool; - typedef std::vector typeVectorVertex; - typedef std::pair typeSimplex; + typedef std::vector typeVectorVertex; + typedef std::pair typeSimplex; std::cout << "********************************************************************" << std::endl; std::cout << "TEST OF RECURSIVE INSERTION" << std::endl; typeST st; @@ -394,7 +398,7 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(NSimplexAndSubfaces_tree_insertion, typeST, list_o // Check it is well inserted BOOST_CHECK(true == returnValue.second); position = 0; - std::sort(SimplexVector1.begin(), SimplexVector1.end(), std::greater()); + std::sort(SimplexVector1.begin(), SimplexVector1.end(), std::greater()); for (auto vertex : st.simplex_vertex_range(returnValue.first)) { // Check returned Simplex_handle std::cout << "vertex = " << vertex << " | vector[" << position << "] = " << SimplexVector1[position] << std::endl; @@ -413,7 +417,7 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(NSimplexAndSubfaces_tree_insertion, typeST, list_o // Check it is well inserted BOOST_CHECK(true == returnValue.second); position = 0; - std::sort(SimplexVector2.begin(), SimplexVector2.end(), std::greater()); + std::sort(SimplexVector2.begin(), SimplexVector2.end(), std::greater()); for (auto vertex : st.simplex_vertex_range(returnValue.first)) { // Check returned Simplex_handle std::cout << "vertex = " << vertex << " | vector[" << position << "] = " << SimplexVector2[position] << std::endl; @@ -432,7 +436,7 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(NSimplexAndSubfaces_tree_insertion, typeST, list_o // Check it is well inserted BOOST_CHECK(true == returnValue.second); position = 0; - std::sort(SimplexVector3.begin(), SimplexVector3.end(), std::greater()); + std::sort(SimplexVector3.begin(), SimplexVector3.end(), std::greater()); for (auto vertex : st.simplex_vertex_range(returnValue.first)) { // Check returned Simplex_handle std::cout << "vertex = " << vertex << " | vector[" << position << "] = " << SimplexVector3[position] << std::endl; @@ -462,7 +466,7 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(NSimplexAndSubfaces_tree_insertion, typeST, list_o // Check it is well inserted BOOST_CHECK(true == returnValue.second); position = 0; - std::sort(SimplexVector5.begin(), SimplexVector5.end(), std::greater()); + std::sort(SimplexVector5.begin(), SimplexVector5.end(), std::greater()); for (auto vertex : st.simplex_vertex_range(returnValue.first)) { // Check returned Simplex_handle std::cout << "vertex = " << vertex << " | vector[" << position << "] = " << SimplexVector5[position] << std::endl; @@ -481,7 +485,7 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(NSimplexAndSubfaces_tree_insertion, typeST, list_o // Check it is well inserted BOOST_CHECK(true == returnValue.second); position = 0; - std::sort(SimplexVector6.begin(), SimplexVector6.end(), std::greater()); + std::sort(SimplexVector6.begin(), SimplexVector6.end(), std::greater()); for (auto vertex : st.simplex_vertex_range(returnValue.first)) { // Check returned Simplex_handle std::cout << "vertex = " << vertex << " | vector[" << position << "] = " << SimplexVector6[position] << std::endl; @@ -504,12 +508,12 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(NSimplexAndSubfaces_tree_insertion, typeST, list_o /* A facet [3,4,5] */ /* A cell [0,1,6,7] */ - typeSimplex simplexPair1 = std::make_pair(SimplexVector1, DEFAULT_FILTRATION_VALUE); - typeSimplex simplexPair2 = std::make_pair(SimplexVector2, DEFAULT_FILTRATION_VALUE); - typeSimplex simplexPair3 = std::make_pair(SimplexVector3, DEFAULT_FILTRATION_VALUE); - typeSimplex simplexPair4 = std::make_pair(SimplexVector4, DEFAULT_FILTRATION_VALUE); - typeSimplex simplexPair5 = std::make_pair(SimplexVector5, DEFAULT_FILTRATION_VALUE); - typeSimplex simplexPair6 = std::make_pair(SimplexVector6, DEFAULT_FILTRATION_VALUE); + typeSimplex simplexPair1 = std::make_pair(SimplexVector1, 0.0); + typeSimplex simplexPair2 = std::make_pair(SimplexVector2, 0.0); + typeSimplex simplexPair3 = std::make_pair(SimplexVector3, 0.0); + typeSimplex simplexPair4 = std::make_pair(SimplexVector4, 0.0); + typeSimplex simplexPair5 = std::make_pair(SimplexVector5, 0.0); + typeSimplex simplexPair6 = std::make_pair(SimplexVector6, 0.0); test_simplex_tree_contains(st, simplexPair1, 6); // (2,1,0) is in position 6 test_simplex_tree_contains(st, simplexPair2, 7); // (3) is in position 7 test_simplex_tree_contains(st, simplexPair3, 8); // (3,0) is in position 8 @@ -600,7 +604,7 @@ void test_cofaces(typeST& st, const std::vector& expected, int di } BOOST_AUTO_TEST_CASE_TEMPLATE(coface_on_simplex_tree, typeST, list_of_tested_variants) { - typedef std::vector typeVectorVertex; + typedef std::vector typeVectorVertex; std::cout << "********************************************************************" << std::endl; std::cout << "TEST COFACE ALGORITHM" << std::endl; typeST st; @@ -629,7 +633,7 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(coface_on_simplex_tree, typeST, list_of_tested_var // FIXME st.set_dimension(3); - std::vector simplex_result; + std::vector simplex_result; std::vector result; std::cout << "First test - Star of (3):" << std::endl; @@ -649,7 +653,7 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(coface_on_simplex_tree, typeST, list_of_tested_var result.push_back(st.find(simplex_result)); simplex_result.clear(); - std::vector vertex = {3}; + std::vector vertex = {3}; test_cofaces(st, vertex, 0, result); vertex.clear(); result.clear(); @@ -699,7 +703,7 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(coface_on_simplex_tree, typeST, list_of_tested_var } BOOST_AUTO_TEST_CASE_TEMPLATE(copy_move_on_simplex_tree, typeST, list_of_tested_variants) { - typedef std::vector typeVectorVertex; + typedef std::vector typeVectorVertex; std::cout << "********************************************************************" << std::endl; std::cout << "TEST COPY MOVE CONSTRUCTORS" << std::endl; typeST st; @@ -771,12 +775,11 @@ void test_simplex_is_vertex(typeST& st, typename typeST::Simplex_handle sh, type BOOST_AUTO_TEST_CASE(non_contiguous) { typedef Simplex_tree<> typeST; - typedef typeST::Vertex_handle Vertex_handle; typedef typeST::Simplex_handle Simplex_handle; std::cout << "********************************************************************" << std::endl; std::cout << "TEST NON-CONTIGUOUS VERTICES" << std::endl; typeST st; - Vertex_handle e[] = {3,-7}; + typeST::Vertex_handle e[] = {3,-7}; std::cout << "Insert" << std::endl; st.insert_simplex_and_subfaces(e); BOOST_CHECK(st.num_vertices() == 2); diff --git a/src/Witness_complex/example/witness_complex_from_file.cpp b/src/Witness_complex/example/witness_complex_from_file.cpp index 53207ad2..392ad55d 100644 --- a/src/Witness_complex/example/witness_complex_from_file.cpp +++ b/src/Witness_complex/example/witness_complex_from_file.cpp @@ -34,7 +34,7 @@ #include #include -typedef std::vector< Vertex_handle > typeVectorVertex; +typedef std::vector< int > typeVectorVertex; typedef std::vector< std::vector > Point_Vector; /** diff --git a/src/Witness_complex/include/gudhi/Witness_complex.h b/src/Witness_complex/include/gudhi/Witness_complex.h index 489cdf11..2cec921a 100644 --- a/src/Witness_complex/include/gudhi/Witness_complex.h +++ b/src/Witness_complex/include/gudhi/Witness_complex.h @@ -72,7 +72,7 @@ class Witness_complex { typedef std::vector< Point_t > Point_Vector; typedef std::vector< Vertex_handle > typeVectorVertex; - typedef std::pair< typeVectorVertex, Filtration_value> typeSimplex; + //typedef std::pair< typeVectorVertex, Filtration_value> typeSimplex; typedef std::pair< Simplex_handle, bool > typePairSimplexBool; typedef int Witness_id; diff --git a/src/Witness_complex/test/simple_witness_complex.cpp b/src/Witness_complex/test/simple_witness_complex.cpp index 03df78ee..adaadfb0 100644 --- a/src/Witness_complex/test/simple_witness_complex.cpp +++ b/src/Witness_complex/test/simple_witness_complex.cpp @@ -33,7 +33,7 @@ #include typedef Gudhi::Simplex_tree<> Simplex_tree; -typedef std::vector< Vertex_handle > typeVectorVertex; +typedef std::vector< int > typeVectorVertex; typedef Gudhi::witness_complex::Witness_complex WitnessComplex; BOOST_AUTO_TEST_CASE(simple_witness_complex) { diff --git a/src/Witness_complex/test/witness_complex_points.cpp b/src/Witness_complex/test/witness_complex_points.cpp index bd3df604..03c9adc0 100644 --- a/src/Witness_complex/test/witness_complex_points.cpp +++ b/src/Witness_complex/test/witness_complex_points.cpp @@ -34,7 +34,7 @@ #include typedef std::vector Point; -typedef std::vector< Vertex_handle > typeVectorVertex; +typedef std::vector< int > typeVectorVertex; typedef Gudhi::Simplex_tree<> Simplex_tree; typedef Gudhi::witness_complex::Witness_complex WitnessComplex; diff --git a/src/common/include/gudhi/reader_utils.h b/src/common/include/gudhi/reader_utils.h index 899f9df6..4154acc9 100644 --- a/src/common/include/gudhi/reader_utils.h +++ b/src/common/include/gudhi/reader_utils.h @@ -77,6 +77,7 @@ inline void read_points(std::string file_name, std::vector< std::vector< double * Every simplex must appear exactly once. * Simplices of dimension more than 1 are ignored. */ +template< typename Graph_t, typename Edge_t, typename Filtration_value, typename Vertex_handle > inline Graph_t read_graph(std::string file_name) { std::ifstream in_(file_name.c_str(), std::ios::in); if (!in_.is_open()) { @@ -130,7 +131,7 @@ inline Graph_t read_graph(std::string file_name) { Graph_t skel_graph(edges.begin(), edges.end(), edges_fil.begin(), vertices.size()); auto vertex_prop = boost::get(vertex_filtration_t(), skel_graph); - boost::graph_traits::vertex_iterator vi, vi_end; + typename boost::graph_traits::vertex_iterator vi, vi_end; auto v_it = vertices.begin(); for (std::tie(vi, vi_end) = boost::vertices(skel_graph); vi != vi_end; ++vi, ++v_it) { boost::put(vertex_prop, *vi, v_it->second); -- cgit v1.2.3 From 64dd2a2b0eec1374ed23ca079f86b312125d03f7 Mon Sep 17 00:00:00 2001 From: vrouvrea Date: Wed, 4 Jan 2017 08:39:18 +0000 Subject: Move bottleneck_chrono in benchmark Add test in cmake for basic example Move CGAL gudhi patches for bottleneck in dedicated directory Fix cpplint syntax issue git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/bottleneck_integration@1920 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: 24671facf791de93dc6fd94bb39ca7362bb22959 --- CMakeLists.txt | 1 + src/Bottleneck_distance/benchmark/CMakeLists.txt | 13 + .../benchmark/bottleneck_chrono.cpp | 62 ++ .../doc/Intro_bottleneck_distance.h | 7 +- src/Bottleneck_distance/example/CMakeLists.txt | 2 + .../example/bottleneck_basic_example.cpp | 26 +- .../example/bottleneck_read_file_example.cpp | 71 ++- src/Bottleneck_distance/include/gudhi/Bottleneck.h | 93 +-- .../include/gudhi/CGAL/Kd_tree.h | 582 ------------------- .../include/gudhi/CGAL/Kd_tree_node.h | 586 ------------------- .../CGAL/Orthogonal_incremental_neighbor_search.h | 621 --------------------- .../include/gudhi/Graph_matching.h | 220 ++++---- .../include/gudhi/Internal_point.h | 69 ++- .../include/gudhi/Neighbors_finder.h | 172 +++--- .../include/gudhi/Persistence_graph.h | 191 ++++--- src/Bottleneck_distance/test/CMakeLists.txt | 1 - src/Bottleneck_distance/test/bottleneck_chrono.cpp | 62 -- .../test/bottleneck_unit_test.cpp | 238 ++++---- src/Tangential_complex/benchmark/CMakeLists.txt | 9 - .../Bottleneck_distance_CGAL_patches.txt | 3 + src/common/include/gudhi_patches/CGAL/Kd_tree.h | 582 +++++++++++++++++++ .../include/gudhi_patches/CGAL/Kd_tree_node.h | 586 +++++++++++++++++++ .../CGAL/Orthogonal_incremental_neighbor_search.h | 621 +++++++++++++++++++++ .../Tangential_complex_CGAL_patches.txt | 82 +++ 24 files changed, 2504 insertions(+), 2396 deletions(-) create mode 100644 src/Bottleneck_distance/benchmark/CMakeLists.txt create mode 100644 src/Bottleneck_distance/benchmark/bottleneck_chrono.cpp delete mode 100644 src/Bottleneck_distance/include/gudhi/CGAL/Kd_tree.h delete mode 100644 src/Bottleneck_distance/include/gudhi/CGAL/Kd_tree_node.h delete mode 100644 src/Bottleneck_distance/include/gudhi/CGAL/Orthogonal_incremental_neighbor_search.h delete mode 100644 src/Bottleneck_distance/test/bottleneck_chrono.cpp create mode 100644 src/common/include/gudhi_patches/Bottleneck_distance_CGAL_patches.txt create mode 100644 src/common/include/gudhi_patches/CGAL/Kd_tree.h create mode 100644 src/common/include/gudhi_patches/CGAL/Kd_tree_node.h create mode 100644 src/common/include/gudhi_patches/CGAL/Orthogonal_incremental_neighbor_search.h create mode 100644 src/common/include/gudhi_patches/Tangential_complex_CGAL_patches.txt (limited to 'CMakeLists.txt') diff --git a/CMakeLists.txt b/CMakeLists.txt index d1a3c7bf..f949a803 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -158,6 +158,7 @@ else() add_subdirectory(src/Tangential_complex/benchmark) add_subdirectory(src/Bottleneck_distance/example) add_subdirectory(src/Bottleneck_distance/test) + add_subdirectory(src/Bottleneck_distance/benchmark) # data points generator add_subdirectory(data/points/generator) diff --git a/src/Bottleneck_distance/benchmark/CMakeLists.txt b/src/Bottleneck_distance/benchmark/CMakeLists.txt new file mode 100644 index 00000000..f70dd8ff --- /dev/null +++ b/src/Bottleneck_distance/benchmark/CMakeLists.txt @@ -0,0 +1,13 @@ +cmake_minimum_required(VERSION 2.6) +project(Bottleneck_distance_benchmark) + + +# requires CGAL 4.8 +# cmake -DCGAL_DIR=~/workspace/CGAL-4.8 ../../.. +if(CGAL_FOUND) + if (NOT CGAL_VERSION VERSION_LESS 4.8.0) + if (EIGEN3_FOUND) + add_executable ( bottleneck_chrono bottleneck_chrono.cpp ) + endif() + endif () +endif() diff --git a/src/Bottleneck_distance/benchmark/bottleneck_chrono.cpp b/src/Bottleneck_distance/benchmark/bottleneck_chrono.cpp new file mode 100644 index 00000000..456c570b --- /dev/null +++ b/src/Bottleneck_distance/benchmark/bottleneck_chrono.cpp @@ -0,0 +1,62 @@ +/* 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: Francois Godi + * + * Copyright (C) 2015 INRIA + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ + +#include +#include +#include +#include + +using namespace Gudhi::persistence_diagram; + + +double upper_bound = 400.; // any real > 0 + +int main() { + std::ofstream result_file; + result_file.open("results.csv", std::ios::out); + + for (int n = 1000; n <= 10000; n += 1000) { + std::uniform_real_distribution unif1(0., upper_bound); + std::uniform_real_distribution unif2(upper_bound / 1000., upper_bound / 100.); + std::default_random_engine re; + std::vector< std::pair > v1, v2; + for (int i = 0; i < n; i++) { + double a = unif1(re); + double b = unif1(re); + double x = unif2(re); + double y = unif2(re); + v1.emplace_back(std::min(a, b), std::max(a, b)); + v2.emplace_back(std::min(a, b) + std::min(x, y), std::max(a, b) + std::max(x, y)); + if (i % 5 == 0) + v1.emplace_back(std::min(a, b), std::min(a, b) + x); + if (i % 3 == 0) + v2.emplace_back(std::max(a, b), std::max(a, b) + y); + } + std::chrono::steady_clock::time_point start = std::chrono::steady_clock::now(); + double b = bottleneck_distance(v1, v2); + std::chrono::steady_clock::time_point end = std::chrono::steady_clock::now(); + typedef std::chrono::duration millisecs_t; + millisecs_t duration(std::chrono::duration_cast(end - start)); + result_file << n << ";" << duration.count() << ";" << b << std::endl; + } + result_file.close(); +} diff --git a/src/Bottleneck_distance/doc/Intro_bottleneck_distance.h b/src/Bottleneck_distance/doc/Intro_bottleneck_distance.h index ebe1123b..21187f9c 100644 --- a/src/Bottleneck_distance/doc/Intro_bottleneck_distance.h +++ b/src/Bottleneck_distance/doc/Intro_bottleneck_distance.h @@ -4,7 +4,7 @@ * * Author: François Godi * - * Copyright (C) 2015 INRIA (France) + * Copyright (C) 2015 INRIA * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by @@ -20,8 +20,8 @@ * along with this program. If not, see . */ -#ifndef DOC_BOTTLENECK_DISTANCE_H_ -#define DOC_BOTTLENECK_DISTANCE_H_ +#ifndef DOC_BOTTLENECK_DISTANCE_INTRO_BOTTLENECK_DISTANCE_H_ +#define DOC_BOTTLENECK_DISTANCE_INTRO_BOTTLENECK_DISTANCE_H_ // needs namespace for Doxygen to link on classes namespace Gudhi { @@ -48,3 +48,4 @@ namespace bottleneck_distance { } // namespace Gudhi +#endif // DOC_BOTTLENECK_DISTANCE_INTRO_BOTTLENECK_DISTANCE_H_ diff --git a/src/Bottleneck_distance/example/CMakeLists.txt b/src/Bottleneck_distance/example/CMakeLists.txt index cd53ccfc..c66623e9 100644 --- a/src/Bottleneck_distance/example/CMakeLists.txt +++ b/src/Bottleneck_distance/example/CMakeLists.txt @@ -8,6 +8,8 @@ if(CGAL_FOUND) if (EIGEN3_FOUND) add_executable (bottleneck_read_file_example bottleneck_read_file_example.cpp) add_executable (bottleneck_basic_example bottleneck_basic_example.cpp) + + add_test(bottleneck_basic_example ${CMAKE_CURRENT_BINARY_DIR}/bottleneck_basic_example) endif() endif () endif() diff --git a/src/Bottleneck_distance/example/bottleneck_basic_example.cpp b/src/Bottleneck_distance/example/bottleneck_basic_example.cpp index 78e00e57..91a7302f 100644 --- a/src/Bottleneck_distance/example/bottleneck_basic_example.cpp +++ b/src/Bottleneck_distance/example/bottleneck_basic_example.cpp @@ -4,7 +4,7 @@ * * Authors: Francois Godi, small modifications by Pawel Dlotko * - * Copyright (C) 2015 INRIA (France) + * Copyright (C) 2015 INRIA * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by @@ -25,24 +25,24 @@ int main() { - std::vector< std::pair > v1, v2; + std::vector< std::pair > v1, v2; - v1.emplace_back(2.7, 3.7); - v1.emplace_back(9.6, 14.); - v1.emplace_back(34.2, 34.974); - v1.emplace_back(3., std::numeric_limits::infinity()); + v1.emplace_back(2.7, 3.7); + v1.emplace_back(9.6, 14.); + v1.emplace_back(34.2, 34.974); + v1.emplace_back(3., std::numeric_limits::infinity()); - v2.emplace_back(2.8, 4.45); - v2.emplace_back(9.5, 14.1); - v2.emplace_back(3.2, std::numeric_limits::infinity()); + v2.emplace_back(2.8, 4.45); + v2.emplace_back(9.5, 14.1); + v2.emplace_back(3.2, std::numeric_limits::infinity()); - double b = Gudhi::persistence_diagram::bottleneck_distance(v1, v2); + double b = Gudhi::persistence_diagram::bottleneck_distance(v1, v2); - std::cout << "Bottleneck distance = " << b << std::endl; + std::cout << "Bottleneck distance = " << b << std::endl; - b = Gudhi::persistence_diagram::bottleneck_distance(v1, v2, 0.1); + b = Gudhi::persistence_diagram::bottleneck_distance(v1, v2, 0.1); - std::cout << "Approx bottleneck distance = " << b << std::endl; + std::cout << "Approx bottleneck distance = " << b << std::endl; } diff --git a/src/Bottleneck_distance/example/bottleneck_read_file_example.cpp b/src/Bottleneck_distance/example/bottleneck_read_file_example.cpp index ceedccc5..4c74b66e 100644 --- a/src/Bottleneck_distance/example/bottleneck_read_file_example.cpp +++ b/src/Bottleneck_distance/example/bottleneck_read_file_example.cpp @@ -4,7 +4,7 @@ * * Authors: Francois Godi, small modifications by Pawel Dlotko * - * Copyright (C) 2015 INRIA (France) + * Copyright (C) 2015 INRIA * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by @@ -28,49 +28,42 @@ #include #include -std::vector< std::pair > read_diagram_from_file( const char* filename ) -{ - std::ifstream in; - in.open( filename ); - std::vector< std::pair > result; - if ( !in.is_open() ) - { - std::cerr << "File : " << filename << " do not exist. The program will now terminate \n"; - throw "File do not exist \n"; - } +std::vector< std::pair > read_diagram_from_file(const char* filename) { + std::ifstream in; + in.open(filename); + std::vector< std::pair > result; + if (!in.is_open()) { + std::cerr << "File : " << filename << " do not exist. The program will now terminate \n"; + throw "File do not exist \n"; + } - std::string line; - while (!in.eof()) - { - getline(in,line); - if ( line.length() != 0 ) - { - std::stringstream lineSS; - lineSS << line; - double beginn, endd; - lineSS >> beginn; - lineSS >> endd; - result.push_back( std::make_pair( beginn , endd ) ); - } + std::string line; + while (!in.eof()) { + getline(in, line); + if (line.length() != 0) { + std::stringstream lineSS; + lineSS << line; + double beginn, endd; + lineSS >> beginn; + lineSS >> endd; + result.push_back(std::make_pair(beginn, endd)); } - in.close(); - return result; -} //read_diagram_from_file + } + in.close(); + return result; +} // read_diagram_from_file -int main( int argc , char** argv ) -{ - if ( argc < 3 ) - { - std::cout << "To run this program please provide as an input two files with persistence diagrams. Each file " << - "should contain a birth-death pair per line. Third, optional parameter is an error bound on a bottleneck" << - " distance (set by default to zero). The program will now terminate \n"; +int main(int argc, char** argv) { + if (argc < 3) { + std::cout << "To run this program please provide as an input two files with persistence diagrams. Each file " << + "should contain a birth-death pair per line. Third, optional parameter is an error bound on a bottleneck" << + " distance (set by default to zero). The program will now terminate \n"; } - std::vector< std::pair< double , double > > diag1 = read_diagram_from_file( argv[1] ); - std::vector< std::pair< double , double > > diag2 = read_diagram_from_file( argv[2] ); + std::vector< std::pair< double, double > > diag1 = read_diagram_from_file(argv[1]); + std::vector< std::pair< double, double > > diag2 = read_diagram_from_file(argv[2]); double tolerance = 0.; - if ( argc == 4 ) - { - tolerance = atof( argv[3] ); + if (argc == 4) { + tolerance = atof(argv[3]); } double b = Gudhi::persistence_diagram::bottleneck_distance(diag1, diag2, tolerance); std::cout << "The distance between the diagrams is : " << b << ". The tolerance is : " << tolerance << std::endl; diff --git a/src/Bottleneck_distance/include/gudhi/Bottleneck.h b/src/Bottleneck_distance/include/gudhi/Bottleneck.h index 42a0d444..2b7e4767 100644 --- a/src/Bottleneck_distance/include/gudhi/Bottleneck.h +++ b/src/Bottleneck_distance/include/gudhi/Bottleneck.h @@ -4,7 +4,7 @@ * * Author: Francois Godi * - * Copyright (C) 2015 INRIA (France) + * Copyright (C) 2015 INRIA * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by @@ -31,62 +31,65 @@ namespace Gudhi { namespace persistence_diagram { double bottleneck_distance_approx(Persistence_graph& g, double e) { - double b_lower_bound = 0.; - double b_upper_bound = g.diameter_bound(); - const double alpha = std::pow(g.size(), 1./5.); - Graph_matching m(g); - Graph_matching biggest_unperfect(g); - while (b_upper_bound - b_lower_bound > 2*e) { - double step = b_lower_bound + (b_upper_bound - b_lower_bound)/alpha; - if(step <= b_lower_bound || step >= b_upper_bound) //Avoid precision problem - break; - m.set_r(step); - while (m.multi_augment()); //compute a maximum matching (in the graph corresponding to the current r) - if (m.perfect()) { - m = biggest_unperfect; - b_upper_bound = step; - } else { - biggest_unperfect = m; - b_lower_bound = step; - } + double b_lower_bound = 0.; + double b_upper_bound = g.diameter_bound(); + const double alpha = std::pow(g.size(), 1. / 5.); + Graph_matching m(g); + Graph_matching biggest_unperfect(g); + while (b_upper_bound - b_lower_bound > 2 * e) { + double step = b_lower_bound + (b_upper_bound - b_lower_bound) / alpha; + if (step <= b_lower_bound || step >= b_upper_bound) // Avoid precision problem + break; + m.set_r(step); + while (m.multi_augment()); // compute a maximum matching (in the graph corresponding to the current r) + if (m.perfect()) { + m = biggest_unperfect; + b_upper_bound = step; + } else { + biggest_unperfect = m; + b_lower_bound = step; } - return (b_lower_bound + b_upper_bound)/2.; + } + return (b_lower_bound + b_upper_bound) / 2.; } double bottleneck_distance_exact(Persistence_graph& g) { - std::vector sd = g.sorted_distances(); - long lower_bound_i = 0; - long upper_bound_i = sd.size()-1; - const double alpha = std::pow(g.size(), 1./5.); - Graph_matching m(g); - Graph_matching biggest_unperfect(g); - while (lower_bound_i != upper_bound_i) { - long step = lower_bound_i + static_cast((upper_bound_i - lower_bound_i - 1)/alpha); - m.set_r(sd.at(step)); - while (m.multi_augment()); //compute a maximum matching (in the graph corresponding to the current r) - if (m.perfect()) { - m = biggest_unperfect; - upper_bound_i = step; - } else { - biggest_unperfect = m; - lower_bound_i = step + 1; - } + std::vector sd = g.sorted_distances(); + long lower_bound_i = 0; + long upper_bound_i = sd.size() - 1; + const double alpha = std::pow(g.size(), 1. / 5.); + Graph_matching m(g); + Graph_matching biggest_unperfect(g); + while (lower_bound_i != upper_bound_i) { + long step = lower_bound_i + static_cast ((upper_bound_i - lower_bound_i - 1) / alpha); + m.set_r(sd.at(step)); + while (m.multi_augment()); // compute a maximum matching (in the graph corresponding to the current r) + if (m.perfect()) { + m = biggest_unperfect; + upper_bound_i = step; + } else { + biggest_unperfect = m; + lower_bound_i = step + 1; } - return sd.at(lower_bound_i); + } + return sd.at(lower_bound_i); } /** \brief Function to use in order to compute the Bottleneck distance between two persistence diagrams (see concepts). - * If the last parameter e is not 0, you get an additive e-approximation, which is a lot faster to compute whatever is e. - * Thus, by default, e is a very small positive double, actually the smallest double possible such that the floating-point inaccuracies don't lead to a failure of the algorithm. + * If the last parameter e is not 0, you get an additive e-approximation, which is a lot faster to compute whatever is + * e. + * Thus, by default, e is a very small positive double, actually the smallest double possible such that the + * floating-point inaccuracies don't lead to a failure of the algorithm. * * \ingroup bottleneck_distance */ template -double bottleneck_distance(const Persistence_diagram1 &diag1, const Persistence_diagram2 &diag2, double e=std::numeric_limits::min()) { - Persistence_graph g(diag1, diag2, e); - if(g.bottleneck_alive() == std::numeric_limits::infinity()) - return std::numeric_limits::infinity(); - return std::max(g.bottleneck_alive(), e == 0. ? bottleneck_distance_exact(g) : bottleneck_distance_approx(g, e)); +double bottleneck_distance(const Persistence_diagram1 &diag1, const Persistence_diagram2 &diag2, + double e = std::numeric_limits::min()) { + Persistence_graph g(diag1, diag2, e); + if (g.bottleneck_alive() == std::numeric_limits::infinity()) + return std::numeric_limits::infinity(); + return std::max(g.bottleneck_alive(), e == 0. ? bottleneck_distance_exact(g) : bottleneck_distance_approx(g, e)); } } // namespace persistence_diagram diff --git a/src/Bottleneck_distance/include/gudhi/CGAL/Kd_tree.h b/src/Bottleneck_distance/include/gudhi/CGAL/Kd_tree.h deleted file mode 100644 index f085b0da..00000000 --- a/src/Bottleneck_distance/include/gudhi/CGAL/Kd_tree.h +++ /dev/null @@ -1,582 +0,0 @@ -// Copyright (c) 2002,2011,2014 Utrecht University (The Netherlands), Max-Planck-Institute Saarbruecken (Germany). -// All rights reserved. -// -// This file is part of CGAL (www.cgal.org). -// You can redistribute it and/or modify it under the terms of the GNU -// General Public License as published by the Free Software Foundation, -// either version 3 of the License, or (at your option) any later version. -// -// Licensees holding a valid commercial license may use this file in -// accordance with the commercial license agreement provided with the software. -// -// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE -// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. -// -// $URL$ -// $Id$ -// -// Author(s) : Hans Tangelder (), -// : Waqar Khan - -#ifndef CGAL_KD_TREE_H -#define CGAL_KD_TREE_H - -#include "Kd_tree_node.h" - -#include -#include -#include - -#include -#include -#include - - -#include -#include -#include - -#ifdef CGAL_HAS_THREADS -#include -#endif - -namespace CGAL { - -//template , class UseExtendedNode = Tag_true > -template , class UseExtendedNode = Tag_true > -class Kd_tree { - -public: - typedef SearchTraits Traits; - typedef Splitter_ Splitter; - typedef typename SearchTraits::Point_d Point_d; - typedef typename Splitter::Container Point_container; - - typedef typename SearchTraits::FT FT; - typedef Kd_tree_node Node; - typedef Kd_tree_leaf_node Leaf_node; - typedef Kd_tree_internal_node Internal_node; - typedef Kd_tree Tree; - typedef Kd_tree Self; - - typedef Node* Node_handle; - typedef const Node* Node_const_handle; - typedef Leaf_node* Leaf_node_handle; - typedef const Leaf_node* Leaf_node_const_handle; - typedef Internal_node* Internal_node_handle; - typedef const Internal_node* Internal_node_const_handle; - typedef typename std::vector::const_iterator Point_d_iterator; - typedef typename std::vector::const_iterator Point_d_const_iterator; - typedef typename Splitter::Separator Separator; - typedef typename std::vector::const_iterator iterator; - typedef typename std::vector::const_iterator const_iterator; - - typedef typename std::vector::size_type size_type; - - typedef typename internal::Get_dimension_tag::Dimension D; - -private: - SearchTraits traits_; - Splitter split; - - - // wokaround for https://svn.boost.org/trac/boost/ticket/9332 -#if (_MSC_VER == 1800) && (BOOST_VERSION == 105500) - std::deque internal_nodes; - std::deque leaf_nodes; -#else - boost::container::deque internal_nodes; - boost::container::deque leaf_nodes; -#endif - - Node_handle tree_root; - - Kd_tree_rectangle* bbox; - std::vector pts; - - // Instead of storing the points in arrays in the Kd_tree_node - // we put all the data in a vector in the Kd_tree. - // and we only store an iterator range in the Kd_tree_node. - // - std::vector data; - - - #ifdef CGAL_HAS_THREADS - mutable CGAL_MUTEX building_mutex;//mutex used to protect const calls inducing build() - #endif - bool built_; - bool removed_; - - // protected copy constructor - Kd_tree(const Tree& tree) - : traits_(tree.traits_),built_(tree.built_) - {}; - - - // Instead of the recursive construction of the tree in the class Kd_tree_node - // we do this in the tree class. The advantage is that we then can optimize - // the allocation of the nodes. - - // The leaf node - Node_handle - create_leaf_node(Point_container& c) - { - Leaf_node node(true , static_cast(c.size())); - std::ptrdiff_t tmp = c.begin() - data.begin(); - node.data = pts.begin() + tmp; - - leaf_nodes.push_back(node); - Leaf_node_handle nh = &leaf_nodes.back(); - - - return nh; - } - - - // The internal node - - Node_handle - create_internal_node(Point_container& c, const Tag_true&) - { - return create_internal_node_use_extension(c); - } - - Node_handle - create_internal_node(Point_container& c, const Tag_false&) - { - return create_internal_node(c); - } - - - - // TODO: Similiar to the leaf_init function above, a part of the code should be - // moved to a the class Kd_tree_node. - // It is not proper yet, but the goal was to see if there is - // a potential performance gain through the Compact_container - Node_handle - create_internal_node_use_extension(Point_container& c) - { - Internal_node node(false); - internal_nodes.push_back(node); - Internal_node_handle nh = &internal_nodes.back(); - - Separator sep; - Point_container c_low(c.dimension(),traits_); - split(sep, c, c_low); - nh->set_separator(sep); - - int cd = nh->cutting_dimension(); - if(!c_low.empty()){ - nh->lower_low_val = c_low.tight_bounding_box().min_coord(cd); - nh->lower_high_val = c_low.tight_bounding_box().max_coord(cd); - } - else{ - nh->lower_low_val = nh->cutting_value(); - nh->lower_high_val = nh->cutting_value(); - } - if(!c.empty()){ - nh->upper_low_val = c.tight_bounding_box().min_coord(cd); - nh->upper_high_val = c.tight_bounding_box().max_coord(cd); - } - else{ - nh->upper_low_val = nh->cutting_value(); - nh->upper_high_val = nh->cutting_value(); - } - - CGAL_assertion(nh->cutting_value() >= nh->lower_low_val); - CGAL_assertion(nh->cutting_value() <= nh->upper_high_val); - - if (c_low.size() > split.bucket_size()){ - nh->lower_ch = create_internal_node_use_extension(c_low); - }else{ - nh->lower_ch = create_leaf_node(c_low); - } - if (c.size() > split.bucket_size()){ - nh->upper_ch = create_internal_node_use_extension(c); - }else{ - nh->upper_ch = create_leaf_node(c); - } - - - - - return nh; - } - - - // Note also that I duplicated the code to get rid if the if's for - // the boolean use_extension which was constant over the construction - Node_handle - create_internal_node(Point_container& c) - { - Internal_node node(false); - internal_nodes.push_back(node); - Internal_node_handle nh = &internal_nodes.back(); - Separator sep; - - Point_container c_low(c.dimension(),traits_); - split(sep, c, c_low); - nh->set_separator(sep); - - if (c_low.size() > split.bucket_size()){ - nh->lower_ch = create_internal_node(c_low); - }else{ - nh->lower_ch = create_leaf_node(c_low); - } - if (c.size() > split.bucket_size()){ - nh->upper_ch = create_internal_node(c); - }else{ - nh->upper_ch = create_leaf_node(c); - } - - - - return nh; - } - - - -public: - - Kd_tree(Splitter s = Splitter(),const SearchTraits traits=SearchTraits()) - : traits_(traits),split(s), built_(false), removed_(false) - {} - - template - Kd_tree(InputIterator first, InputIterator beyond, - Splitter s = Splitter(),const SearchTraits traits=SearchTraits()) - : traits_(traits),split(s), built_(false), removed_(false) - { - pts.insert(pts.end(), first, beyond); - } - - bool empty() const { - return pts.empty(); - } - - void - build() - { - // This function is not ready to be called when a tree already exists, one - // must call invalidate_built() first. - CGAL_assertion(!is_built()); - CGAL_assertion(!removed_); - const Point_d& p = *pts.begin(); - typename SearchTraits::Construct_cartesian_const_iterator_d ccci=traits_.construct_cartesian_const_iterator_d_object(); - int dim = static_cast(std::distance(ccci(p), ccci(p,0))); - - data.reserve(pts.size()); - for(unsigned int i = 0; i < pts.size(); i++){ - data.push_back(&pts[i]); - } - Point_container c(dim, data.begin(), data.end(),traits_); - bbox = new Kd_tree_rectangle(c.bounding_box()); - if (c.size() <= split.bucket_size()){ - tree_root = create_leaf_node(c); - }else { - tree_root = create_internal_node(c, UseExtendedNode()); - } - - //Reorder vector for spatial locality - std::vector ptstmp; - ptstmp.resize(pts.size()); - for (std::size_t i = 0; i < pts.size(); ++i){ - ptstmp[i] = *data[i]; - } - for(std::size_t i = 0; i < leaf_nodes.size(); ++i){ - std::ptrdiff_t tmp = leaf_nodes[i].begin() - pts.begin(); - leaf_nodes[i].data = ptstmp.begin() + tmp; - } - pts.swap(ptstmp); - - data.clear(); - - built_ = true; - } - -private: - //any call to this function is for the moment not threadsafe - void const_build() const { - #ifdef CGAL_HAS_THREADS - //this ensure that build() will be called once - CGAL_SCOPED_LOCK(building_mutex); - if(!is_built()) - #endif - const_cast(this)->build(); //THIS IS NOT THREADSAFE - } -public: - - bool is_built() const - { - return built_; - } - - void invalidate_built() - { - if(removed_){ - // Walk the tree to collect the remaining points. - // Writing directly to pts would likely work, but better be safe. - std::vector ptstmp; - //ptstmp.resize(root()->num_items()); - root()->tree_items(std::back_inserter(ptstmp)); - pts.swap(ptstmp); - removed_=false; - CGAL_assertion(is_built()); // the rest of the cleanup must happen - } - if(is_built()){ - internal_nodes.clear(); - leaf_nodes.clear(); - data.clear(); - delete bbox; - built_ = false; - } - } - - void clear() - { - invalidate_built(); - pts.clear(); - removed_ = false; - } - - void - insert(const Point_d& p) - { - invalidate_built(); - pts.push_back(p); - } - - template - void - insert(InputIterator first, InputIterator beyond) - { - invalidate_built(); - pts.insert(pts.end(),first, beyond); - } - -private: - struct Equal_by_coordinates { - SearchTraits const* traits; - Point_d const* pp; - bool operator()(Point_d const&q) const { - typename SearchTraits::Construct_cartesian_const_iterator_d ccci=traits->construct_cartesian_const_iterator_d_object(); - return std::equal(ccci(*pp), ccci(*pp,0), ccci(q)); - } - }; - Equal_by_coordinates equal_by_coordinates(Point_d const&p){ - Equal_by_coordinates ret = { &traits(), &p }; - return ret; - } - -public: - void - remove(const Point_d& p) - { - remove(p, equal_by_coordinates(p)); - } - - template - void - remove(const Point_d& p, Equal const& equal_to_p) - { -#if 0 - // This code could have quadratic runtime. - if (!is_built()) { - std::vector::iterator pi = std::find(pts.begin(), pts.end(), p); - // Precondition: the point must be there. - CGAL_assertion (pi != pts.end()); - pts.erase(pi); - return; - } -#endif - bool success = remove_(p, 0, false, 0, false, root(), equal_to_p); - CGAL_assertion(success); - - // Do not set the flag is the tree has been cleared. - if(is_built()) - removed_ |= success; - } -private: - template - bool remove_(const Point_d& p, - Internal_node_handle grandparent, bool parent_islower, - Internal_node_handle parent, bool islower, - Node_handle node, Equal const& equal_to_p) { - // Recurse to locate the point - if (!node->is_leaf()) { - Internal_node_handle newparent = static_cast(node); - // FIXME: This should be if(xcutting_dimension()] <= newparent->cutting_value()) { - if (remove_(p, parent, islower, newparent, true, newparent->lower(), equal_to_p)) - return true; - } - //if (traits().construct_cartesian_const_iterator_d_object()(p)[newparent->cutting_dimension()] >= newparent->cutting_value()) - return remove_(p, parent, islower, newparent, false, newparent->upper(), equal_to_p); - - CGAL_assertion(false); // Point was not found - } - - // Actual removal - Leaf_node_handle lnode = static_cast(node); - if (lnode->size() > 1) { - iterator pi = std::find_if(lnode->begin(), lnode->end(), equal_to_p); - // FIXME: we should ensure this never happens - if (pi == lnode->end()) return false; - iterator lasti = lnode->end() - 1; - if (pi != lasti) { - // Hack to get a non-const iterator - std::iter_swap(pts.begin()+(pi-pts.begin()), pts.begin()+(lasti-pts.begin())); - } - lnode->drop_last_point(); - } else if (!equal_to_p(*lnode->begin())) { - // FIXME: we should ensure this never happens - return false; - } else if (grandparent) { - Node_handle brother = islower ? parent->upper() : parent->lower(); - if (parent_islower) - grandparent->set_lower(brother); - else - grandparent->set_upper(brother); - } else if (parent) { - tree_root = islower ? parent->upper() : parent->lower(); - } else { - clear(); - } - return true; - } - -public: - //For efficiency; reserve the size of the points vectors in advance (if the number of points is already known). - void reserve(size_t size) - { - pts.reserve(size); - } - - //Get the capacity of the underlying points vector. - size_t capacity() - { - return pts.capacity(); - } - - - template - OutputIterator - search(OutputIterator it, const FuzzyQueryItem& q) const - { - if(! pts.empty()){ - - if(! is_built()){ - const_build(); - } - Kd_tree_rectangle b(*bbox); - return tree_root->search(it,q,b); - } - return it; - } - - - template - boost::optional - search_any_point(const FuzzyQueryItem& q) const - { - if(! pts.empty()){ - - if(! is_built()){ - const_build(); - } - Kd_tree_rectangle b(*bbox); - return tree_root->search_any_point(q,b); - } - return boost::none; - } - - - ~Kd_tree() { - if(is_built()){ - delete bbox; - } - } - - - const SearchTraits& - traits() const - { - return traits_; - } - - Node_const_handle - root() const - { - if(! is_built()){ - const_build(); - } - return tree_root; - } - - Node_handle - root() - { - if(! is_built()){ - build(); - } - return tree_root; - } - - void - print() const - { - if(! is_built()){ - const_build(); - } - root()->print(); - } - - const Kd_tree_rectangle& - bounding_box() const - { - if(! is_built()){ - const_build(); - } - return *bbox; - } - - const_iterator - begin() const - { - return pts.begin(); - } - - const_iterator - end() const - { - return pts.end(); - } - - size_type - size() const - { - return pts.size(); - } - - // Print statistics of the tree. - std::ostream& - statistics(std::ostream& s) const - { - if(! is_built()){ - const_build(); - } - s << "Tree statistics:" << std::endl; - s << "Number of items stored: " - << root()->num_items() << std::endl; - s << "Number of nodes: " - << root()->num_nodes() << std::endl; - s << " Tree depth: " << root()->depth() << std::endl; - return s; - } - - -}; - -} // namespace CGAL - -#endif // CGAL_KD_TREE_H diff --git a/src/Bottleneck_distance/include/gudhi/CGAL/Kd_tree_node.h b/src/Bottleneck_distance/include/gudhi/CGAL/Kd_tree_node.h deleted file mode 100644 index 909ee260..00000000 --- a/src/Bottleneck_distance/include/gudhi/CGAL/Kd_tree_node.h +++ /dev/null @@ -1,586 +0,0 @@ -// Copyright (c) 2002,2011 Utrecht University (The Netherlands). -// All rights reserved. -// -// This file is part of CGAL (www.cgal.org). -// You can redistribute it and/or modify it under the terms of the GNU -// General Public License as published by the Free Software Foundation, -// either version 3 of the License, or (at your option) any later version. -// -// Licensees holding a valid commercial license may use this file in -// accordance with the commercial license agreement provided with the software. -// -// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE -// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. -// -// $URL$ -// $Id$ -// -// -// Authors : Hans Tangelder () - -#ifndef CGAL_KD_TREE_NODE_H -#define CGAL_KD_TREE_NODE_H - -#include "CGAL/Splitters.h" - -#include -#include - -namespace CGAL { - - template - class Kd_tree; - - template < class TreeTraits, class Splitter, class UseExtendedNode > - class Kd_tree_node { - - friend class Kd_tree; - - typedef typename Kd_tree::Node_handle Node_handle; - typedef typename Kd_tree::Node_const_handle Node_const_handle; - typedef typename Kd_tree::Internal_node_handle Internal_node_handle; - typedef typename Kd_tree::Internal_node_const_handle Internal_node_const_handle; - typedef typename Kd_tree::Leaf_node_handle Leaf_node_handle; - typedef typename Kd_tree::Leaf_node_const_handle Leaf_node_const_handle; - typedef typename TreeTraits::Point_d Point_d; - - typedef typename TreeTraits::FT FT; - typedef typename Kd_tree::Separator Separator; - typedef typename Kd_tree::Point_d_iterator Point_d_iterator; - typedef typename Kd_tree::iterator iterator; - typedef typename Kd_tree::D D; - - bool leaf; - - public : - Kd_tree_node(bool leaf_) - :leaf(leaf_){} - - bool is_leaf() const{ - return leaf; - } - - std::size_t - num_items() const - { - if (is_leaf()){ - Leaf_node_const_handle node = - static_cast(this); - return node->size(); - } - else { - Internal_node_const_handle node = - static_cast(this); - return node->lower()->num_items() + node->upper()->num_items(); - } - } - - std::size_t - num_nodes() const - { - if (is_leaf()) return 1; - else { - Internal_node_const_handle node = - static_cast(this); - return node->lower()->num_nodes() + node->upper()->num_nodes(); - } - } - - int - depth(const int current_max_depth) const - { - if (is_leaf()){ - return current_max_depth; - } - else { - Internal_node_const_handle node = - static_cast(this); - return - (std::max)( node->lower()->depth(current_max_depth + 1), - node->upper()->depth(current_max_depth + 1)); - } - } - - int - depth() const - { - return depth(1); - } - - template - OutputIterator - tree_items(OutputIterator it) const { - if (is_leaf()) { - Leaf_node_const_handle node = - static_cast(this); - if (node->size()>0) - for (iterator i=node->begin(); i != node->end(); i++) - {*it=*i; ++it;} - } - else { - Internal_node_const_handle node = - static_cast(this); - it=node->lower()->tree_items(it); - it=node->upper()->tree_items(it); - } - return it; - } - - - boost::optional - any_tree_item() const { - boost::optional result = boost::none; - if (is_leaf()) { - Leaf_node_const_handle node = - static_cast(this); - if (node->size()>0){ - return boost::make_optional(*(node->begin())); - } - } - else { - Internal_node_const_handle node = - static_cast(this); - result = node->lower()->any_tree_item(); - if(! result){ - result = node->upper()->any_tree_item(); - } - } - return result; - } - - - void - indent(int d) const - { - for(int i = 0; i < d; i++){ - std::cout << " "; - } - } - - - void - print(int d = 0) const - { - if (is_leaf()) { - Leaf_node_const_handle node = - static_cast(this); - indent(d); - std::cout << "leaf" << std::endl; - if (node->size()>0) - for (iterator i=node->begin(); i != node->end(); i++) - {indent(d);std::cout << *i << std::endl;} - } - else { - Internal_node_const_handle node = - static_cast(this); - indent(d); - std::cout << "lower tree" << std::endl; - node->lower()->print(d+1); - indent(d); - std::cout << "separator: dim = " << node->cutting_dimension() << " val = " << node->cutting_value() << std::endl; - indent(d); - std::cout << "upper tree" << std::endl; - node->upper()->print(d+1); - } - } - - - template - OutputIterator - search(OutputIterator it, const FuzzyQueryItem& q, - Kd_tree_rectangle& b) const - { - if (is_leaf()) { - Leaf_node_const_handle node = - static_cast(this); - if (node->size()>0) - for (iterator i=node->begin(); i != node->end(); i++) - if (q.contains(*i)) - {*it++=*i;} - } - else { - Internal_node_const_handle node = - static_cast(this); - // after splitting b denotes the lower part of b - Kd_tree_rectangle b_upper(b); - b.split(b_upper, node->cutting_dimension(), - node->cutting_value()); - - if (q.outer_range_contains(b)) - it=node->lower()->tree_items(it); - else - if (q.inner_range_intersects(b)) - it=node->lower()->search(it,q,b); - if (q.outer_range_contains(b_upper)) - it=node->upper()->tree_items(it); - else - if (q.inner_range_intersects(b_upper)) - it=node->upper()->search(it,q,b_upper); - }; - return it; - } - - - template - boost::optional - search_any_point(const FuzzyQueryItem& q, - Kd_tree_rectangle& b) const - { - boost::optional result = boost::none; - if (is_leaf()) { - Leaf_node_const_handle node = - static_cast(this); - if (node->size()>0) - for (iterator i=node->begin(); i != node->end(); i++) - if (q.contains(*i)) - { result = *i; break; } - } - else { - Internal_node_const_handle node = - static_cast(this); - // after splitting b denotes the lower part of b - Kd_tree_rectangle b_upper(b); - b.split(b_upper, node->cutting_dimension(), - node->cutting_value()); - - if (q.outer_range_contains(b)){ - result = node->lower()->any_tree_item(); - }else{ - if (q.inner_range_intersects(b)){ - result = node->lower()->search_any_point(q,b); - } - } - if(result){ - return result; - } - if (q.outer_range_contains(b_upper)){ - result = node->upper()->any_tree_item(); - }else{ - if (q.inner_range_intersects(b_upper)) - result = node->upper()->search_any_point(q,b_upper); - } - } - return result; - } - - }; - - - template < class TreeTraits, class Splitter, class UseExtendedNode > - class Kd_tree_leaf_node : public Kd_tree_node< TreeTraits, Splitter, UseExtendedNode >{ - - friend class Kd_tree; - - typedef typename Kd_tree::iterator iterator; - typedef Kd_tree_node< TreeTraits, Splitter, UseExtendedNode> Base; - typedef typename TreeTraits::Point_d Point_d; - - private: - - // private variables for leaf nodes - boost::int32_t n; // denotes number of items in a leaf node - iterator data; // iterator to data in leaf node - - - public: - - // default constructor - Kd_tree_leaf_node() - {} - - Kd_tree_leaf_node(bool leaf_ ) - : Base(leaf_) - {} - - Kd_tree_leaf_node(bool leaf_,unsigned int n_ ) - : Base(leaf_), n(n_) - {} - - // members for all nodes - - // members for leaf nodes only - inline - unsigned int - size() const - { - return n; - } - - inline - iterator - begin() const - { - return data; - } - - inline - iterator - end() const - { - return data + n; - } - - inline - void - drop_last_point() - { - --n; - } - - }; //leaf node - - - - template < class TreeTraits, class Splitter, class UseExtendedNode> - class Kd_tree_internal_node : public Kd_tree_node< TreeTraits, Splitter, UseExtendedNode >{ - - friend class Kd_tree; - - typedef Kd_tree_node< TreeTraits, Splitter, UseExtendedNode> Base; - typedef typename Kd_tree::Node_handle Node_handle; - typedef typename Kd_tree::Node_const_handle Node_const_handle; - - typedef typename TreeTraits::FT FT; - typedef typename Kd_tree::Separator Separator; - - private: - - // private variables for internal nodes - boost::int32_t cut_dim; - FT cut_val; - Node_handle lower_ch, upper_ch; - - - // private variables for extended internal nodes - FT upper_low_val; - FT upper_high_val; - FT lower_low_val; - FT lower_high_val; - - - public: - - // default constructor - Kd_tree_internal_node() - {} - - Kd_tree_internal_node(bool leaf_) - : Base(leaf_) - {} - - - // members for internal node and extended internal node - - inline - Node_const_handle - lower() const - { - return lower_ch; - } - - inline - Node_const_handle - upper() const - { - return upper_ch; - } - - inline - Node_handle - lower() - { - return lower_ch; - } - - inline - Node_handle - upper() - { - return upper_ch; - } - - inline - void - set_lower(Node_handle nh) - { - lower_ch = nh; - } - - inline - void - set_upper(Node_handle nh) - { - upper_ch = nh; - } - - // inline Separator& separator() {return sep; } - // use instead - inline - void set_separator(Separator& sep){ - cut_dim = sep.cutting_dimension(); - cut_val = sep.cutting_value(); - } - - inline - FT - cutting_value() const - { - return cut_val; - } - - inline - int - cutting_dimension() const - { - return cut_dim; - } - - // members for extended internal node only - inline - FT - upper_low_value() const - { - return upper_low_val; - } - - inline - FT - upper_high_value() const - { - return upper_high_val; - } - - inline - FT - lower_low_value() const - { - return lower_low_val; - } - - inline - FT - lower_high_value() const - { - return lower_high_val; - } - - /*Separator& - separator() - { - return Separator(cutting_dimension,cutting_value); - }*/ - - - };//internal node - - template < class TreeTraits, class Splitter> - class Kd_tree_internal_node : public Kd_tree_node< TreeTraits, Splitter, Tag_false >{ - - friend class Kd_tree; - - typedef Kd_tree_node< TreeTraits, Splitter, Tag_false> Base; - typedef typename Kd_tree::Node_handle Node_handle; - typedef typename Kd_tree::Node_const_handle Node_const_handle; - - typedef typename TreeTraits::FT FT; - typedef typename Kd_tree::Separator Separator; - - private: - - // private variables for internal nodes - boost::uint8_t cut_dim; - FT cut_val; - - Node_handle lower_ch, upper_ch; - - public: - - // default constructor - Kd_tree_internal_node() - {} - - Kd_tree_internal_node(bool leaf_) - : Base(leaf_) - {} - - - // members for internal node and extended internal node - - inline - Node_const_handle - lower() const - { - return lower_ch; - } - - inline - Node_const_handle - upper() const - { - return upper_ch; - } - - inline - Node_handle - lower() - { - return lower_ch; - } - - inline - Node_handle - upper() - { - return upper_ch; - } - - inline - void - set_lower(Node_handle nh) - { - lower_ch = nh; - } - - inline - void - set_upper(Node_handle nh) - { - upper_ch = nh; - } - - // inline Separator& separator() {return sep; } - // use instead - - inline - void set_separator(Separator& sep){ - cut_dim = sep.cutting_dimension(); - cut_val = sep.cutting_value(); - } - - inline - FT - cutting_value() const - { - return cut_val; - } - - inline - int - cutting_dimension() const - { - return cut_dim; - } - - /* Separator& - separator() - { - return Separator(cutting_dimension,cutting_value); - }*/ - - - };//internal node - - - -} // namespace CGAL -#endif // CGAL_KDTREE_NODE_H diff --git a/src/Bottleneck_distance/include/gudhi/CGAL/Orthogonal_incremental_neighbor_search.h b/src/Bottleneck_distance/include/gudhi/CGAL/Orthogonal_incremental_neighbor_search.h deleted file mode 100644 index dbe707ed..00000000 --- a/src/Bottleneck_distance/include/gudhi/CGAL/Orthogonal_incremental_neighbor_search.h +++ /dev/null @@ -1,621 +0,0 @@ -// Copyright (c) 2002,2011 Utrecht University (The Netherlands). -// All rights reserved. -// -// This file is part of CGAL (www.cgal.org). -// You can redistribute it and/or modify it under the terms of the GNU -// General Public License as published by the Free Software Foundation, -// either version 3 of the License, or (at your option) any later version. -// -// Licensees holding a valid commercial license may use this file in -// accordance with the commercial license agreement provided with the software. -// -// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE -// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. -// -// $URL$ -// $Id$ -// -// -// Author(s) : Hans Tangelder () - -#ifndef CGAL_ORTHOGONAL_INCREMENTAL_NEIGHBOR_SEARCH -#define CGAL_ORTHOGONAL_INCREMENTAL_NEIGHBOR_SEARCH - -#include "CGAL/Kd_tree.h" - -#include -#include -#include -#include -#include -#include - -namespace CGAL { - - template ::type, - class Splitter_ = Sliding_midpoint, - class Tree_= Kd_tree > - class Orthogonal_incremental_neighbor_search { - - public: - typedef Splitter_ Splitter; - typedef Tree_ Tree; - typedef Distance_ Distance; - typedef typename SearchTraits::Point_d Point_d; - typedef typename Distance::Query_item Query_item; - typedef typename SearchTraits::FT FT; - typedef typename Tree::Point_d_iterator Point_d_iterator; - typedef typename Tree::Node_const_handle Node_const_handle; - - typedef std::pair Point_with_transformed_distance; - typedef CGAL::cpp11::tuple > Node_with_distance; - typedef std::vector Node_with_distance_vector; - typedef std::vector Point_with_transformed_distance_vector; - - template - struct Object_wrapper - { - T object; - Object_wrapper(const T& t):object(t){} - const T& operator* () const { return object; } - const T* operator-> () const { return &object; } - }; - - class Iterator_implementation { - SearchTraits traits; - public: - - int number_of_neighbours_computed; - int number_of_internal_nodes_visited; - int number_of_leaf_nodes_visited; - int number_of_items_visited; - - private: - - typedef std::vector Distance_vector; - - Distance_vector dists; - - Distance Orthogonal_distance_instance; - - FT multiplication_factor; - - Query_item query_point; - - FT distance_to_root; - - bool search_nearest_neighbour; - - FT rd; - - - class Priority_higher { - public: - - bool search_nearest; - - Priority_higher(bool search_the_nearest_neighbour) - : search_nearest(search_the_nearest_neighbour) - {} - - //highest priority is smallest distance - bool - operator() (Node_with_distance* n1, Node_with_distance* n2) const - { - return (search_nearest) ? (CGAL::cpp11::get<1>(*n1) > CGAL::cpp11::get<1>(*n2)) : (CGAL::cpp11::get<1>(*n2) > CGAL::cpp11::get<1>(*n1)); - } - }; - - class Distance_smaller { - - public: - - bool search_nearest; - - Distance_smaller(bool search_the_nearest_neighbour) - : search_nearest(search_the_nearest_neighbour) - {} - - //highest priority is smallest distance - bool operator() (Point_with_transformed_distance* p1, Point_with_transformed_distance* p2) const - { - return (search_nearest) ? (p1->second > p2->second) : (p2->second > p1->second); - } - }; - - - std::priority_queue PriorityQueue; - - public: - std::priority_queue Item_PriorityQueue; - - - public: - - int reference_count; - - - - // constructor - Iterator_implementation(const Tree& tree,const Query_item& q, const Distance& tr, - FT Eps=FT(0.0), bool search_nearest=true) - : traits(tree.traits()),number_of_neighbours_computed(0), number_of_internal_nodes_visited(0), - number_of_leaf_nodes_visited(0), number_of_items_visited(0), - Orthogonal_distance_instance(tr), multiplication_factor(Orthogonal_distance_instance.transformed_distance(FT(1.0)+Eps)), - query_point(q), search_nearest_neighbour(search_nearest), - PriorityQueue(Priority_higher(search_nearest)), Item_PriorityQueue(Distance_smaller(search_nearest)), - reference_count(1) - - - { - if (tree.empty()) return; - - typename SearchTraits::Construct_cartesian_const_iterator_d ccci=traits.construct_cartesian_const_iterator_d_object(); - int dim = static_cast(std::distance(ccci(q), ccci(q,0))); - - dists.resize(dim); - for(int i=0 ; i(*The_Root); - Compute_the_next_nearest_neighbour(); - } - else{ - distance_to_root= - Orthogonal_distance_instance.max_distance_to_rectangle(q, - tree.bounding_box(), dists); - Node_with_distance *The_Root = new Node_with_distance(tree.root(), - distance_to_root, dists); - PriorityQueue.push(The_Root); - - // rd is the distance of the top of the priority queue to q - rd=CGAL::cpp11::get<1>(*The_Root); - Compute_the_next_furthest_neighbour(); - } - - - } - - // * operator - const Point_with_transformed_distance& - operator* () const - { - return *(Item_PriorityQueue.top()); - } - - // prefix operator - Iterator_implementation& - operator++() - { - Delete_the_current_item_top(); - if(search_nearest_neighbour) - Compute_the_next_nearest_neighbour(); - else - Compute_the_next_furthest_neighbour(); - return *this; - } - - // postfix operator - Object_wrapper - operator++(int) - { - Object_wrapper result( *(Item_PriorityQueue.top()) ); - ++*this; - return result; - } - - // Print statistics of the general priority search process. - std::ostream& - statistics (std::ostream& s) const { - s << "Orthogonal priority search statistics:" - << std::endl; - s << "Number of internal nodes visited:" - << number_of_internal_nodes_visited << std::endl; - s << "Number of leaf nodes visited:" - << number_of_leaf_nodes_visited << std::endl; - s << "Number of items visited:" - << number_of_items_visited << std::endl; - s << "Number of neighbours computed:" - << number_of_neighbours_computed << std::endl; - return s; - } - - - //destructor - ~Iterator_implementation() - { - while (!PriorityQueue.empty()) { - Node_with_distance* The_top=PriorityQueue.top(); - PriorityQueue.pop(); - delete The_top; - } - while (!Item_PriorityQueue.empty()) { - Point_with_transformed_distance* The_top=Item_PriorityQueue.top(); - Item_PriorityQueue.pop(); - delete The_top; - } - } - - private: - - void - Delete_the_current_item_top() - { - Point_with_transformed_distance* The_item_top=Item_PriorityQueue.top(); - Item_PriorityQueue.pop(); - delete The_item_top; - } - - void - Compute_the_next_nearest_neighbour() - { - // compute the next item - bool next_neighbour_found=false; - if (!(Item_PriorityQueue.empty())) { - next_neighbour_found= - (multiplication_factor*rd > Item_PriorityQueue.top()->second); - } - typename SearchTraits::Construct_cartesian_const_iterator_d construct_it=traits.construct_cartesian_const_iterator_d_object(); - typename SearchTraits::Cartesian_const_iterator_d query_point_it = construct_it(query_point); - // otherwise browse the tree further - while ((!next_neighbour_found) && (!PriorityQueue.empty())) { - Node_with_distance* The_node_top=PriorityQueue.top(); - Node_const_handle N= CGAL::cpp11::get<0>(*The_node_top); - dists = CGAL::cpp11::get<2>(*The_node_top); - PriorityQueue.pop(); - delete The_node_top; - FT copy_rd=rd; - while (!(N->is_leaf())) { // compute new distance - typename Tree::Internal_node_const_handle node = - static_cast(N); - number_of_internal_nodes_visited++; - int new_cut_dim=node->cutting_dimension(); - FT new_rd,dst = dists[new_cut_dim]; - FT val = *(query_point_it + new_cut_dim); - FT diff1 = val - node->upper_low_value(); - FT diff2 = val - node->lower_high_value(); - if (diff1 + diff2 < FT(0.0)) { - new_rd= - Orthogonal_distance_instance.new_distance(copy_rd,dst,diff1,new_cut_dim); - - CGAL_assertion(new_rd >= copy_rd); - dists[new_cut_dim] = diff1; - Node_with_distance *Upper_Child = - new Node_with_distance(node->upper(), new_rd, dists); - PriorityQueue.push(Upper_Child); - dists[new_cut_dim] = dst; - N=node->lower(); - - } - else { // compute new distance - new_rd=Orthogonal_distance_instance.new_distance(copy_rd,dst,diff2,new_cut_dim); - CGAL_assertion(new_rd >= copy_rd); - dists[new_cut_dim] = diff2; - Node_with_distance *Lower_Child = - new Node_with_distance(node->lower(), new_rd, dists); - PriorityQueue.push(Lower_Child); - dists[new_cut_dim] = dst; - N=node->upper(); - } - } - // n is a leaf - typename Tree::Leaf_node_const_handle node = - static_cast(N); - number_of_leaf_nodes_visited++; - if (node->size() > 0) { - for (typename Tree::iterator it=node->begin(); it != node->end(); it++) { - number_of_items_visited++; - FT distance_to_query_point= - Orthogonal_distance_instance.transformed_distance(query_point,*it); - Point_with_transformed_distance *NN_Candidate= - new Point_with_transformed_distance(*it,distance_to_query_point); - Item_PriorityQueue.push(NN_Candidate); - } - // old top of PriorityQueue has been processed, - // hence update rd - - if (!(PriorityQueue.empty())) { - rd = CGAL::cpp11::get<1>(*PriorityQueue.top()); - next_neighbour_found = - (multiplication_factor*rd > - Item_PriorityQueue.top()->second); - } - else // priority queue empty => last neighbour found - { - next_neighbour_found=true; - } - - number_of_neighbours_computed++; - } - } // next_neighbour_found or priority queue is empty - // in the latter case also the item priority quee is empty - } - - - void - Compute_the_next_furthest_neighbour() - { - // compute the next item - bool next_neighbour_found=false; - if (!(Item_PriorityQueue.empty())) { - next_neighbour_found= - (rd < multiplication_factor*Item_PriorityQueue.top()->second); - } - typename SearchTraits::Construct_cartesian_const_iterator_d construct_it=traits.construct_cartesian_const_iterator_d_object(); - typename SearchTraits::Cartesian_const_iterator_d query_point_it = construct_it(query_point); - // otherwise browse the tree further - while ((!next_neighbour_found) && (!PriorityQueue.empty())) { - Node_with_distance* The_node_top=PriorityQueue.top(); - Node_const_handle N= CGAL::cpp11::get<0>(*The_node_top); - dists = CGAL::cpp11::get<2>(*The_node_top); - PriorityQueue.pop(); - delete The_node_top; - FT copy_rd=rd; - while (!(N->is_leaf())) { // compute new distance - typename Tree::Internal_node_const_handle node = - static_cast(N); - number_of_internal_nodes_visited++; - int new_cut_dim=node->cutting_dimension(); - FT new_rd,dst = dists[new_cut_dim]; - FT val = *(query_point_it + new_cut_dim); - FT diff1 = val - node->upper_low_value(); - FT diff2 = val - node->lower_high_value(); - if (diff1 + diff2 < FT(0.0)) { - diff1 = val - node->upper_high_value(); - new_rd= - Orthogonal_distance_instance.new_distance(copy_rd,dst,diff1,new_cut_dim); - Node_with_distance *Lower_Child = - new Node_with_distance(node->lower(), copy_rd, dists); - PriorityQueue.push(Lower_Child); - N=node->upper(); - dists[new_cut_dim] = diff1; - copy_rd=new_rd; - - } - else { // compute new distance - diff2 = val - node->lower_low_value(); - new_rd=Orthogonal_distance_instance.new_distance(copy_rd,dst,diff2,new_cut_dim); - Node_with_distance *Upper_Child = - new Node_with_distance(node->upper(), copy_rd, dists); - PriorityQueue.push(Upper_Child); - N=node->lower(); - dists[new_cut_dim] = diff2; - copy_rd=new_rd; - } - } - // n is a leaf - typename Tree::Leaf_node_const_handle node = - static_cast(N); - number_of_leaf_nodes_visited++; - if (node->size() > 0) { - for (typename Tree::iterator it=node->begin(); it != node->end(); it++) { - number_of_items_visited++; - FT distance_to_query_point= - Orthogonal_distance_instance.transformed_distance(query_point,*it); - Point_with_transformed_distance *NN_Candidate= - new Point_with_transformed_distance(*it,distance_to_query_point); - Item_PriorityQueue.push(NN_Candidate); - } - // old top of PriorityQueue has been processed, - // hence update rd - - if (!(PriorityQueue.empty())) { - rd = CGAL::cpp11::get<1>(*PriorityQueue.top()); - next_neighbour_found = - (multiplication_factor*rd < - Item_PriorityQueue.top()->second); - } - else // priority queue empty => last neighbour found - { - next_neighbour_found=true; - } - - number_of_neighbours_computed++; - } - } // next_neighbour_found or priority queue is empty - // in the latter case also the item priority quee is empty - } - }; // class Iterator_implementaion - - - - - - - - - - public: - class iterator; - typedef iterator const_iterator; - - // constructor - Orthogonal_incremental_neighbor_search(const Tree& tree, - const Query_item& q, FT Eps = FT(0.0), - bool search_nearest=true, const Distance& tr=Distance()) - : m_tree(tree),m_query(q),m_dist(tr),m_Eps(Eps),m_search_nearest(search_nearest) - {} - - iterator - begin() const - { - return iterator(m_tree,m_query,m_dist,m_Eps,m_search_nearest); - } - - iterator - end() const - { - return iterator(); - } - - std::ostream& - statistics(std::ostream& s) - { - begin()->statistics(s); - return s; - } - - - - - class iterator { - - public: - - typedef std::input_iterator_tag iterator_category; - typedef Point_with_transformed_distance value_type; - typedef Point_with_transformed_distance* pointer; - typedef const Point_with_transformed_distance& reference; - typedef std::size_t size_type; - typedef std::ptrdiff_t difference_type; - typedef int distance_type; - - //class Iterator_implementation; - Iterator_implementation *Ptr_implementation; - - - public: - - // default constructor - iterator() - : Ptr_implementation(0) - {} - - int - the_number_of_items_visited() - { - return Ptr_implementation->number_of_items_visited; - } - - // constructor - iterator(const Tree& tree,const Query_item& q, const Distance& tr=Distance(), FT eps=FT(0.0), - bool search_nearest=true) - : Ptr_implementation(new Iterator_implementation(tree, q, tr, eps, search_nearest)) - {} - - // copy constructor - iterator(const iterator& Iter) - { - Ptr_implementation = Iter.Ptr_implementation; - if (Ptr_implementation != 0) Ptr_implementation->reference_count++; - } - - iterator& operator=(const iterator& Iter) - { - if (Ptr_implementation != Iter.Ptr_implementation){ - if (Ptr_implementation != 0 && --(Ptr_implementation->reference_count)==0) { - delete Ptr_implementation; - } - Ptr_implementation = Iter.Ptr_implementation; - if (Ptr_implementation != 0) Ptr_implementation->reference_count++; - } - return *this; - } - - - const Point_with_transformed_distance& - operator* () const - { - return *(*Ptr_implementation); - } - - // -> operator - const Point_with_transformed_distance* - operator-> () const - { - return &*(*Ptr_implementation); - } - - // prefix operator - iterator& - operator++() - { - ++(*Ptr_implementation); - return *this; - } - - // postfix operator - Object_wrapper - operator++(int) - { - return (*Ptr_implementation)++; - } - - - bool - operator==(const iterator& It) const - { - if ( - ((Ptr_implementation == 0) || - Ptr_implementation->Item_PriorityQueue.empty()) && - ((It.Ptr_implementation == 0) || - It.Ptr_implementation->Item_PriorityQueue.empty()) - ) - return true; - // else - return (Ptr_implementation == It.Ptr_implementation); - } - - bool - operator!=(const iterator& It) const - { - return !(*this == It); - } - - std::ostream& - statistics (std::ostream& s) - { - Ptr_implementation->statistics(s); - return s; - } - - ~iterator() - { - if (Ptr_implementation != 0) { - Ptr_implementation->reference_count--; - if (Ptr_implementation->reference_count==0) { - delete Ptr_implementation; - Ptr_implementation = 0; - } - } - } - - - }; // class iterator - - //data members - const Tree& m_tree; - Query_item m_query; - Distance m_dist; - FT m_Eps; - bool m_search_nearest; - }; // class - - template - void swap (typename Orthogonal_incremental_neighbor_search::iterator& x, - typename Orthogonal_incremental_neighbor_search::iterator& y) - { - typename Orthogonal_incremental_neighbor_search::iterator::Iterator_implementation - *tmp = x.Ptr_implementation; - x.Ptr_implementation = y.Ptr_implementation; - y.Ptr_implementation = tmp; - } - -} // namespace CGAL - -#endif // CGAL_ORTHOGONAL_INCREMENTAL_NEIGHBOR_SEARCH_H diff --git a/src/Bottleneck_distance/include/gudhi/Graph_matching.h b/src/Bottleneck_distance/include/gudhi/Graph_matching.h index e9f455d7..253c89b4 100644 --- a/src/Bottleneck_distance/include/gudhi/Graph_matching.h +++ b/src/Bottleneck_distance/include/gudhi/Graph_matching.h @@ -4,7 +4,7 @@ * * Author: Francois Godi * - * Copyright (C) 2015 INRIA (France) + * Copyright (C) 2015 INRIA * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by @@ -34,141 +34,141 @@ namespace persistence_diagram { * \ingroup bottleneck_distance */ class Graph_matching { -public: - /** \internal \brief Constructor constructing an empty matching. */ - explicit Graph_matching(Persistence_graph &g); - /** \internal \brief Copy operator. */ - Graph_matching& operator=(const Graph_matching& m); - /** \internal \brief Is the matching perfect ? */ - bool perfect() const; - /** \internal \brief Augments the matching with a maximal set of edge-disjoint shortest augmenting paths. */ - bool multi_augment(); - /** \internal \brief Sets the maximum length of the edges allowed to be added in the matching, 0 initially. */ - void set_r(double r); - -private: - Persistence_graph& g; - double r; - /** \internal \brief Given a point from V, provides its matched point in U, null_point_index() if there isn't. */ - std::vector v_to_u; - /** \internal \brief All the unmatched points in U. */ - std::list unmatched_in_u; - - /** \internal \brief Provides a Layered_neighbors_finder dividing the graph in layers. Basically a BFS. */ - Layered_neighbors_finder layering() const; - /** \internal \brief Augments the matching with a simple path no longer than max_depth. Basically a DFS. */ - bool augment(Layered_neighbors_finder & layered_nf, int u_start_index, int max_depth); - /** \internal \brief Update the matching with the simple augmenting path given as parameter. */ - void update(std::vector & path); + public: + /** \internal \brief Constructor constructing an empty matching. */ + explicit Graph_matching(Persistence_graph &g); + /** \internal \brief Copy operator. */ + Graph_matching& operator=(const Graph_matching& m); + /** \internal \brief Is the matching perfect ? */ + bool perfect() const; + /** \internal \brief Augments the matching with a maximal set of edge-disjoint shortest augmenting paths. */ + bool multi_augment(); + /** \internal \brief Sets the maximum length of the edges allowed to be added in the matching, 0 initially. */ + void set_r(double r); + + private: + Persistence_graph& g; + double r; + /** \internal \brief Given a point from V, provides its matched point in U, null_point_index() if there isn't. */ + std::vector v_to_u; + /** \internal \brief All the unmatched points in U. */ + std::list unmatched_in_u; + + /** \internal \brief Provides a Layered_neighbors_finder dividing the graph in layers. Basically a BFS. */ + Layered_neighbors_finder layering() const; + /** \internal \brief Augments the matching with a simple path no longer than max_depth. Basically a DFS. */ + bool augment(Layered_neighbors_finder & layered_nf, int u_start_index, int max_depth); + /** \internal \brief Update the matching with the simple augmenting path given as parameter. */ + void update(std::vector & path); }; inline Graph_matching::Graph_matching(Persistence_graph& g) : g(g), r(0.), v_to_u(g.size(), null_point_index()), unmatched_in_u() { - for (int u_point_index = 0; u_point_index < g.size(); ++u_point_index) - unmatched_in_u.emplace_back(u_point_index); + for (int u_point_index = 0; u_point_index < g.size(); ++u_point_index) + unmatched_in_u.emplace_back(u_point_index); } inline Graph_matching& Graph_matching::operator=(const Graph_matching& m) { - g = m.g; - r = m.r; - v_to_u = m.v_to_u; - unmatched_in_u = m.unmatched_in_u; - return *this; + g = m.g; + r = m.r; + v_to_u = m.v_to_u; + unmatched_in_u = m.unmatched_in_u; + return *this; } inline bool Graph_matching::perfect() const { - return unmatched_in_u.empty(); + return unmatched_in_u.empty(); } inline bool Graph_matching::multi_augment() { - if (perfect()) - return false; - Layered_neighbors_finder layered_nf(layering()); - int max_depth = layered_nf.vlayers_number()*2 - 1; - double rn = sqrt(g.size()); - // verification of a necessary criterion in order to shortcut if possible - if (max_depth <0 || (unmatched_in_u.size() > rn && max_depth >= rn)) - return false; - bool successful = false; - std::list tries(unmatched_in_u); - for (auto it = tries.cbegin(); it != tries.cend(); it++) - // 'augment' has side-effects which have to be always executed, don't change order - successful = augment(layered_nf, *it, max_depth) || successful; - return successful; + if (perfect()) + return false; + Layered_neighbors_finder layered_nf(layering()); + int max_depth = layered_nf.vlayers_number()*2 - 1; + double rn = sqrt(g.size()); + // verification of a necessary criterion in order to shortcut if possible + if (max_depth < 0 || (unmatched_in_u.size() > rn && max_depth >= rn)) + return false; + bool successful = false; + std::list tries(unmatched_in_u); + for (auto it = tries.cbegin(); it != tries.cend(); it++) + // 'augment' has side-effects which have to be always executed, don't change order + successful = augment(layered_nf, *it, max_depth) || successful; + return successful; } inline void Graph_matching::set_r(double r) { - this->r = r; + this->r = r; } inline bool Graph_matching::augment(Layered_neighbors_finder & layered_nf, int u_start_index, int max_depth) { - //V vertices have at most one successor, thus when we backtrack from U we can directly pop_back 2 vertices. - std::vector path; - path.emplace_back(u_start_index); - do { - if (static_cast(path.size()) > max_depth) { - path.pop_back(); - path.pop_back(); - } - if (path.empty()) - return false; - path.emplace_back(layered_nf.pull_near(path.back(), static_cast(path.size())/2)); - while (path.back() == null_point_index()) { - path.pop_back(); - path.pop_back(); - if (path.empty()) - return false; - path.pop_back(); - path.emplace_back(layered_nf.pull_near(path.back(), path.size() / 2)); - } - path.emplace_back(v_to_u.at(path.back())); - } while (path.back() != null_point_index()); - //if v_to_u.at(path.back()) has no successor, path.back() is an exposed vertex - path.pop_back(); - update(path); - return true; + // V vertices have at most one successor, thus when we backtrack from U we can directly pop_back 2 vertices. + std::vector path; + path.emplace_back(u_start_index); + do { + if (static_cast (path.size()) > max_depth) { + path.pop_back(); + path.pop_back(); + } + if (path.empty()) + return false; + path.emplace_back(layered_nf.pull_near(path.back(), static_cast (path.size()) / 2)); + while (path.back() == null_point_index()) { + path.pop_back(); + path.pop_back(); + if (path.empty()) + return false; + path.pop_back(); + path.emplace_back(layered_nf.pull_near(path.back(), path.size() / 2)); + } + path.emplace_back(v_to_u.at(path.back())); + } while (path.back() != null_point_index()); + // if v_to_u.at(path.back()) has no successor, path.back() is an exposed vertex + path.pop_back(); + update(path); + return true; } inline Layered_neighbors_finder Graph_matching::layering() const { - std::list u_vertices(unmatched_in_u); - std::list v_vertices; - Neighbors_finder nf(g, r); - for (int v_point_index = 0; v_point_index < g.size(); ++v_point_index) - nf.add(v_point_index); - Layered_neighbors_finder layered_nf(g, r); - for(int layer = 0; !u_vertices.empty(); layer++) { - // one layer is one step in the BFS - for (auto it1 = u_vertices.cbegin(); it1 != u_vertices.cend(); ++it1) { - std::vector u_succ(nf.pull_all_near(*it1)); - for (auto it2 = u_succ.begin(); it2 != u_succ.end(); ++it2) { - layered_nf.add(*it2, layer); - v_vertices.emplace_back(*it2); - } - } - // When the above for finishes, we have progress of one half-step (from U to V) in the BFS - u_vertices.clear(); - bool end = false; - for (auto it = v_vertices.cbegin(); it != v_vertices.cend(); it++) - if (v_to_u.at(*it) == null_point_index()) - // we stop when a nearest exposed V vertex (from U exposed vertices) has been found - end = true; - else - u_vertices.emplace_back(v_to_u.at(*it)); - // When the above for finishes, we have progress of one half-step (from V to U) in the BFS - if (end) - return layered_nf; - v_vertices.clear(); + std::list u_vertices(unmatched_in_u); + std::list v_vertices; + Neighbors_finder nf(g, r); + for (int v_point_index = 0; v_point_index < g.size(); ++v_point_index) + nf.add(v_point_index); + Layered_neighbors_finder layered_nf(g, r); + for (int layer = 0; !u_vertices.empty(); layer++) { + // one layer is one step in the BFS + for (auto it1 = u_vertices.cbegin(); it1 != u_vertices.cend(); ++it1) { + std::vector u_succ(nf.pull_all_near(*it1)); + for (auto it2 = u_succ.begin(); it2 != u_succ.end(); ++it2) { + layered_nf.add(*it2, layer); + v_vertices.emplace_back(*it2); + } } - return layered_nf; + // When the above for finishes, we have progress of one half-step (from U to V) in the BFS + u_vertices.clear(); + bool end = false; + for (auto it = v_vertices.cbegin(); it != v_vertices.cend(); it++) + if (v_to_u.at(*it) == null_point_index()) + // we stop when a nearest exposed V vertex (from U exposed vertices) has been found + end = true; + else + u_vertices.emplace_back(v_to_u.at(*it)); + // When the above for finishes, we have progress of one half-step (from V to U) in the BFS + if (end) + return layered_nf; + v_vertices.clear(); + } + return layered_nf; } inline void Graph_matching::update(std::vector& path) { - unmatched_in_u.remove(path.front()); - for (auto it = path.cbegin(); it != path.cend(); ++it) { - // Be careful, the iterator is incremented twice each time - int tmp = *it; - v_to_u[*(++it)] = tmp; - } + unmatched_in_u.remove(path.front()); + for (auto it = path.cbegin(); it != path.cend(); ++it) { + // Be careful, the iterator is incremented twice each time + int tmp = *it; + v_to_u[*(++it)] = tmp; + } } diff --git a/src/Bottleneck_distance/include/gudhi/Internal_point.h b/src/Bottleneck_distance/include/gudhi/Internal_point.h index 70342d0c..0b2d26fe 100644 --- a/src/Bottleneck_distance/include/gudhi/Internal_point.h +++ b/src/Bottleneck_distance/include/gudhi/Internal_point.h @@ -4,7 +4,7 @@ * * Author: Francois Godi * - * Copyright (C) 2015 INRIA (France) + * Copyright (C) 2015 INRIA * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by @@ -32,39 +32,60 @@ int null_point_index(); /** \internal \typedef \brief Internal_point is the internal points representation, indexes used outside. */ struct Internal_point { - double vec[2]; - int point_index; - Internal_point() {} - Internal_point(double x, double y, int p_i) { vec[0]=x; vec[1]=y; point_index = p_i; } - double x() const { return vec[ 0 ]; } - double y() const { return vec[ 1 ]; } - double& x() { return vec[ 0 ]; } - double& y() { return vec[ 1 ]; } - bool operator==(const Internal_point& p) const - { - return point_index==p.point_index; - } - bool operator!=(const Internal_point& p) const { return !(*this == p); } + double vec[2]; + int point_index; + + Internal_point() { } + + Internal_point(double x, double y, int p_i) { + vec[0] = x; + vec[1] = y; + point_index = p_i; + } + + double x() const { + return vec[ 0 ]; + } + + double y() const { + return vec[ 1 ]; + } + + double& x() { + return vec[ 0 ]; + } + + double& y() { + return vec[ 1 ]; + } + + bool operator==(const Internal_point& p) const { + return point_index == p.point_index; + } + + bool operator!=(const Internal_point& p) const { + return !(*this == p); + } }; inline int null_point_index() { - return -1; + return -1; } struct Construct_coord_iterator { - typedef const double* result_type; - const double* operator()(const Internal_point& p) const - { return p.vec; } - const double* operator()(const Internal_point& p, int) const - { return p.vec+2; } + typedef const double* result_type; + + const double* operator()(const Internal_point& p) const { + return p.vec; + } + + const double* operator()(const Internal_point& p, int) const { + return p.vec + 2; + } }; } // namespace persistence_diagram } // namespace Gudhi - - - - #endif // INTERNAL_POINT_H_ diff --git a/src/Bottleneck_distance/include/gudhi/Neighbors_finder.h b/src/Bottleneck_distance/include/gudhi/Neighbors_finder.h index 792925b7..96ece360 100644 --- a/src/Bottleneck_distance/include/gudhi/Neighbors_finder.h +++ b/src/Bottleneck_distance/include/gudhi/Neighbors_finder.h @@ -4,7 +4,7 @@ * * Author: Francois Godi * - * Copyright (C) 2015 INRIA (France) + * Copyright (C) 2015 INRIA * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by @@ -24,10 +24,9 @@ #define NEIGHBORS_FINDER_H_ // Inclusion order is important for CGAL patch -#include "CGAL/Kd_tree_node.h" -#include "CGAL/Kd_tree.h" -#include "CGAL/Orthogonal_k_neighbor_search.h" - +#include +#include +#include #include #include @@ -43,123 +42,126 @@ namespace persistence_diagram { /** \internal \brief data structure used to find any point (including projections) in V near to a query point from U * (which can be a projection). * - * V points have to be added manually using their index and before the first pull. A neighbor pulled is automatically removed. + * V points have to be added manually using their index and before the first pull. A neighbor pulled is automatically + * removed. * * \ingroup bottleneck_distance */ class Neighbors_finder { - - typedef CGAL::Dimension_tag<2> D; - typedef CGAL::Search_traits Traits; - typedef CGAL::Weighted_Minkowski_distance Distance; - typedef CGAL::Orthogonal_k_neighbor_search K_neighbor_search; - typedef K_neighbor_search::Tree Kd_tree; - -public: - /** \internal \brief Constructor taking the near distance definition as parameter. */ - Neighbors_finder(const Persistence_graph& g, double r); - /** \internal \brief A point added will be possibly pulled. */ - void add(int v_point_index); - /** \internal \brief Returns and remove a V point near to the U point given as parameter, null_point_index() if there isn't such a point. */ - int pull_near(int u_point_index); - /** \internal \brief Returns and remove all the V points near to the U point given as parameter. */ - std::vector pull_all_near(int u_point_index); - -private: - const Persistence_graph& g; - const double r; - Kd_tree kd_t; - std::unordered_set projections_f; + typedef CGAL::Dimension_tag<2> D; + typedef CGAL::Search_traits Traits; + typedef CGAL::Weighted_Minkowski_distance Distance; + typedef CGAL::Orthogonal_k_neighbor_search K_neighbor_search; + typedef K_neighbor_search::Tree Kd_tree; + + public: + /** \internal \brief Constructor taking the near distance definition as parameter. */ + Neighbors_finder(const Persistence_graph& g, double r); + /** \internal \brief A point added will be possibly pulled. */ + void add(int v_point_index); + /** \internal \brief Returns and remove a V point near to the U point given as parameter, null_point_index() if + * there isn't such a point. */ + int pull_near(int u_point_index); + /** \internal \brief Returns and remove all the V points near to the U point given as parameter. */ + std::vector pull_all_near(int u_point_index); + + private: + const Persistence_graph& g; + const double r; + Kd_tree kd_t; + std::unordered_set projections_f; }; /** \internal \brief data structure used to find any point (including projections) in V near to a query point from U * (which can be a projection) in a layered graph layer given as parmeter. * - * V points have to be added manually using their index and before the first pull. A neighbor pulled is automatically removed. + * V points have to be added manually using their index and before the first pull. A neighbor pulled is automatically + * removed. * * \ingroup bottleneck_distance */ class Layered_neighbors_finder { -public: - /** \internal \brief Constructor taking the near distance definition as parameter. */ - Layered_neighbors_finder(const Persistence_graph& g, double r); - /** \internal \brief A point added will be possibly pulled. */ - void add(int v_point_index, int vlayer); - /** \internal \brief Returns and remove a V point near to the U point given as parameter, null_point_index() if there isn't such a point. */ - int pull_near(int u_point_index, int vlayer); - /** \internal \brief Returns the number of layers. */ - int vlayers_number() const; - -private: - const Persistence_graph& g; - const double r; - std::vector> neighbors_finder; + public: + /** \internal \brief Constructor taking the near distance definition as parameter. */ + Layered_neighbors_finder(const Persistence_graph& g, double r); + /** \internal \brief A point added will be possibly pulled. */ + void add(int v_point_index, int vlayer); + /** \internal \brief Returns and remove a V point near to the U point given as parameter, null_point_index() if + * there isn't such a point. */ + int pull_near(int u_point_index, int vlayer); + /** \internal \brief Returns the number of layers. */ + int vlayers_number() const; + + private: + const Persistence_graph& g; + const double r; + std::vector> neighbors_finder; }; inline Neighbors_finder::Neighbors_finder(const Persistence_graph& g, double r) : g(g), r(r), kd_t(), projections_f() { } inline void Neighbors_finder::add(int v_point_index) { - if (g.on_the_v_diagonal(v_point_index)) - projections_f.emplace(v_point_index); - else - kd_t.insert(g.get_v_point(v_point_index)); + if (g.on_the_v_diagonal(v_point_index)) + projections_f.emplace(v_point_index); + else + kd_t.insert(g.get_v_point(v_point_index)); } inline int Neighbors_finder::pull_near(int u_point_index) { - int tmp; - int c = g.corresponding_point_in_v(u_point_index); - if (g.on_the_u_diagonal(u_point_index) && !projections_f.empty()){ - //Any pair of projection is at distance 0 - tmp = *projections_f.cbegin(); - projections_f.erase(tmp); - } - else if (projections_f.count(c) && (g.distance(u_point_index, c) <= r)){ - //Is the query point near to its projection ? - tmp = c; - projections_f.erase(tmp); - } - else{ - //Is the query point near to a V point in the plane ? - Internal_point u_point = g.get_u_point(u_point_index); - std::array w = { {1., 1.} }; - K_neighbor_search search(kd_t, u_point, 1, 0., true, Distance(0, 2, w.begin(), w.end())); - auto it = search.begin(); - if(it==search.end() || g.distance(u_point_index, it->first.point_index) > r) - return null_point_index(); - tmp = it->first.point_index; - kd_t.remove(g.get_v_point(tmp)); - } - return tmp; + int tmp; + int c = g.corresponding_point_in_v(u_point_index); + if (g.on_the_u_diagonal(u_point_index) && !projections_f.empty()) { + // Any pair of projection is at distance 0 + tmp = *projections_f.cbegin(); + projections_f.erase(tmp); + } else if (projections_f.count(c) && (g.distance(u_point_index, c) <= r)) { + // Is the query point near to its projection ? + tmp = c; + projections_f.erase(tmp); + } else { + // Is the query point near to a V point in the plane ? + Internal_point u_point = g.get_u_point(u_point_index); + std::array w = { + {1., 1.} + }; + K_neighbor_search search(kd_t, u_point, 1, 0., true, Distance(0, 2, w.begin(), w.end())); + auto it = search.begin(); + if (it == search.end() || g.distance(u_point_index, it->first.point_index) > r) + return null_point_index(); + tmp = it->first.point_index; + kd_t.remove(g.get_v_point(tmp)); + } + return tmp; } inline std::vector Neighbors_finder::pull_all_near(int u_point_index) { - std::vector all_pull; - int last_pull = pull_near(u_point_index); - while (last_pull != null_point_index()) { - all_pull.emplace_back(last_pull); - last_pull = pull_near(u_point_index); - } - return all_pull; + std::vector all_pull; + int last_pull = pull_near(u_point_index); + while (last_pull != null_point_index()) { + all_pull.emplace_back(last_pull); + last_pull = pull_near(u_point_index); + } + return all_pull; } inline Layered_neighbors_finder::Layered_neighbors_finder(const Persistence_graph& g, double r) : g(g), r(r), neighbors_finder() { } inline void Layered_neighbors_finder::add(int v_point_index, int vlayer) { - for (int l = neighbors_finder.size(); l <= vlayer; l++) - neighbors_finder.emplace_back(std::unique_ptr(new Neighbors_finder(g, r))); - neighbors_finder.at(vlayer)->add(v_point_index); + for (int l = neighbors_finder.size(); l <= vlayer; l++) + neighbors_finder.emplace_back(std::unique_ptr(new Neighbors_finder(g, r))); + neighbors_finder.at(vlayer)->add(v_point_index); } inline int Layered_neighbors_finder::pull_near(int u_point_index, int vlayer) { - if (static_cast (neighbors_finder.size()) <= vlayer) - return null_point_index(); - return neighbors_finder.at(vlayer)->pull_near(u_point_index); + if (static_cast (neighbors_finder.size()) <= vlayer) + return null_point_index(); + return neighbors_finder.at(vlayer)->pull_near(u_point_index); } inline int Layered_neighbors_finder::vlayers_number() const { - return static_cast(neighbors_finder.size()); + return static_cast (neighbors_finder.size()); } } // namespace persistence_diagram diff --git a/src/Bottleneck_distance/include/gudhi/Persistence_graph.h b/src/Bottleneck_distance/include/gudhi/Persistence_graph.h index 45a4d586..3a4a5fec 100644 --- a/src/Bottleneck_distance/include/gudhi/Persistence_graph.h +++ b/src/Bottleneck_distance/include/gudhi/Persistence_graph.h @@ -4,7 +4,7 @@ * * Author: Francois Godi * - * Copyright (C) 2015 INRIA (France) + * Copyright (C) 2015 INRIA * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by @@ -31,147 +31,144 @@ namespace Gudhi { namespace persistence_diagram { - /** \internal \brief Structure representing an euclidean bipartite graph containing * the points from the two persistence diagrams (including the projections). * * \ingroup bottleneck_distance */ class Persistence_graph { -public: - /** \internal \brief Constructor taking 2 Persistence_Diagrams (concept) as parameters. */ - template - Persistence_graph(const Persistence_diagram1& diag1, const Persistence_diagram2& diag2, double e); - /** \internal \brief Is the given point from U the projection of a point in V ? */ - bool on_the_u_diagonal(int u_point_index) const; - /** \internal \brief Is the given point from V the projection of a point in U ? */ - bool on_the_v_diagonal(int v_point_index) const; - /** \internal \brief Given a point from V, returns the corresponding (projection or projector) point in U. */ - int corresponding_point_in_u(int v_point_index) const; - /** \internal \brief Given a point from U, returns the corresponding (projection or projector) point in V. */ - int corresponding_point_in_v(int u_point_index) const; - /** \internal \brief Given a point from U and a point from V, returns the distance between those points. */ - double distance(int u_point_index, int v_point_index) const; - /** \internal \brief Returns size = |U| = |V|. */ - int size() const; - /** \internal \brief Is there as many infinite points (alive components) in both diagrams ? */ - double bottleneck_alive() const; - /** \internal \brief Returns the O(n^2) sorted distances between the points. */ - std::vector sorted_distances() const; - /** \internal \brief Returns an upper bound for the diameter of the convex hull of all non infinite points */ - double diameter_bound() const; - /** \internal \brief Returns the corresponding internal point */ - Internal_point get_u_point(int u_point_index) const; - /** \internal \brief Returns the corresponding internal point */ - Internal_point get_v_point(int v_point_index) const; - -private: - std::vector u; - std::vector v; - double b_alive; + public: + /** \internal \brief Constructor taking 2 Persistence_Diagrams (concept) as parameters. */ + template + Persistence_graph(const Persistence_diagram1& diag1, const Persistence_diagram2& diag2, double e); + /** \internal \brief Is the given point from U the projection of a point in V ? */ + bool on_the_u_diagonal(int u_point_index) const; + /** \internal \brief Is the given point from V the projection of a point in U ? */ + bool on_the_v_diagonal(int v_point_index) const; + /** \internal \brief Given a point from V, returns the corresponding (projection or projector) point in U. */ + int corresponding_point_in_u(int v_point_index) const; + /** \internal \brief Given a point from U, returns the corresponding (projection or projector) point in V. */ + int corresponding_point_in_v(int u_point_index) const; + /** \internal \brief Given a point from U and a point from V, returns the distance between those points. */ + double distance(int u_point_index, int v_point_index) const; + /** \internal \brief Returns size = |U| = |V|. */ + int size() const; + /** \internal \brief Is there as many infinite points (alive components) in both diagrams ? */ + double bottleneck_alive() const; + /** \internal \brief Returns the O(n^2) sorted distances between the points. */ + std::vector sorted_distances() const; + /** \internal \brief Returns an upper bound for the diameter of the convex hull of all non infinite points */ + double diameter_bound() const; + /** \internal \brief Returns the corresponding internal point */ + Internal_point get_u_point(int u_point_index) const; + /** \internal \brief Returns the corresponding internal point */ + Internal_point get_v_point(int v_point_index) const; + + private: + std::vector u; + std::vector v; + double b_alive; }; template Persistence_graph::Persistence_graph(const Persistence_diagram1 &diag1, const Persistence_diagram2 &diag2, double e) - : u(), v(), b_alive(0.) -{ - std::vector u_alive; - std::vector v_alive; - for (auto it = std::begin(diag1); it != std::end(diag1); ++it){ - if(std::get<1>(*it) == std::numeric_limits::infinity()) - u_alive.push_back(std::get<0>(*it)); - else if (std::get<1>(*it) - std::get<0>(*it) > e) - u.push_back(Internal_point(std::get<0>(*it), std::get<1>(*it), u.size())); - } - for (auto it = std::begin(diag2); it != std::end(diag2); ++it){ - if(std::get<1>(*it) == std::numeric_limits::infinity()) - v_alive.push_back(std::get<0>(*it)); - else if (std::get<1>(*it) - std::get<0>(*it) > e) - v.push_back(Internal_point(std::get<0>(*it), std::get<1>(*it), v.size())); - } - if (u.size() < v.size()) - swap(u, v); - std::sort(u_alive.begin(), u_alive.end()); - std::sort(v_alive.begin(), v_alive.end()); - if(u_alive.size() != v_alive.size()) - b_alive = std::numeric_limits::infinity(); - else for(auto it_u=u_alive.cbegin(), it_v=v_alive.cbegin(); it_u != u_alive.cend(); ++it_u, ++it_v) - b_alive = std::max(b_alive, std::fabs(*it_u - *it_v)); + : u(), v(), b_alive(0.) { + std::vector u_alive; + std::vector v_alive; + for (auto it = std::begin(diag1); it != std::end(diag1); ++it) { + if (std::get<1>(*it) == std::numeric_limits::infinity()) + u_alive.push_back(std::get<0>(*it)); + else if (std::get<1>(*it) - std::get<0>(*it) > e) + u.push_back(Internal_point(std::get<0>(*it), std::get<1>(*it), u.size())); + } + for (auto it = std::begin(diag2); it != std::end(diag2); ++it) { + if (std::get<1>(*it) == std::numeric_limits::infinity()) + v_alive.push_back(std::get<0>(*it)); + else if (std::get<1>(*it) - std::get<0>(*it) > e) + v.push_back(Internal_point(std::get<0>(*it), std::get<1>(*it), v.size())); + } + if (u.size() < v.size()) + swap(u, v); + std::sort(u_alive.begin(), u_alive.end()); + std::sort(v_alive.begin(), v_alive.end()); + if (u_alive.size() != v_alive.size()) + b_alive = std::numeric_limits::infinity(); + else for (auto it_u = u_alive.cbegin(), it_v = v_alive.cbegin(); it_u != u_alive.cend(); ++it_u, ++it_v) + b_alive = std::max(b_alive, std::fabs(*it_u - *it_v)); } inline bool Persistence_graph::on_the_u_diagonal(int u_point_index) const { - return u_point_index >= static_cast (u.size()); + return u_point_index >= static_cast (u.size()); } inline bool Persistence_graph::on_the_v_diagonal(int v_point_index) const { - return v_point_index >= static_cast (v.size()); + return v_point_index >= static_cast (v.size()); } inline int Persistence_graph::corresponding_point_in_u(int v_point_index) const { - return on_the_v_diagonal(v_point_index) ? - v_point_index - static_cast (v.size()) : v_point_index + static_cast (u.size()); + return on_the_v_diagonal(v_point_index) ? + v_point_index - static_cast (v.size()) : v_point_index + static_cast (u.size()); } inline int Persistence_graph::corresponding_point_in_v(int u_point_index) const { - return on_the_u_diagonal(u_point_index) ? - u_point_index - static_cast (u.size()) : u_point_index + static_cast (v.size()); + return on_the_u_diagonal(u_point_index) ? + u_point_index - static_cast (u.size()) : u_point_index + static_cast (v.size()); } inline double Persistence_graph::distance(int u_point_index, int v_point_index) const { - if (on_the_u_diagonal(u_point_index) && on_the_v_diagonal(v_point_index)) - return 0.; - Internal_point p_u = get_u_point(u_point_index); - Internal_point p_v = get_v_point(v_point_index); - return std::max(std::fabs(p_u.x() - p_v.x()), std::fabs(p_u.y() - p_v.y())); + if (on_the_u_diagonal(u_point_index) && on_the_v_diagonal(v_point_index)) + return 0.; + Internal_point p_u = get_u_point(u_point_index); + Internal_point p_v = get_v_point(v_point_index); + return std::max(std::fabs(p_u.x() - p_v.x()), std::fabs(p_u.y() - p_v.y())); } inline int Persistence_graph::size() const { - return static_cast (u.size() + v.size()); + return static_cast (u.size() + v.size()); } -inline double Persistence_graph::bottleneck_alive() const{ - return b_alive; +inline double Persistence_graph::bottleneck_alive() const { + return b_alive; } inline std::vector Persistence_graph::sorted_distances() const { - std::vector distances; - distances.push_back(0.); //for empty diagrams - for (int u_point_index = 0; u_point_index < size(); ++u_point_index){ - distances.push_back(distance(u_point_index, corresponding_point_in_v(u_point_index))); - for (int v_point_index = 0; v_point_index < size(); ++v_point_index) - distances.push_back(distance(u_point_index, v_point_index)); - } - std::sort(distances.begin(), distances.end()); - return distances; + std::vector distances; + distances.push_back(0.); // for empty diagrams + for (int u_point_index = 0; u_point_index < size(); ++u_point_index) { + distances.push_back(distance(u_point_index, corresponding_point_in_v(u_point_index))); + for (int v_point_index = 0; v_point_index < size(); ++v_point_index) + distances.push_back(distance(u_point_index, v_point_index)); + } + std::sort(distances.begin(), distances.end()); + return distances; } inline Internal_point Persistence_graph::get_u_point(int u_point_index) const { - if (!on_the_u_diagonal(u_point_index)) - return u.at(u_point_index); - Internal_point projector = v.at(corresponding_point_in_v(u_point_index)); - double m = (projector.x() + projector.y()) / 2.; - return Internal_point(m,m,u_point_index); + if (!on_the_u_diagonal(u_point_index)) + return u.at(u_point_index); + Internal_point projector = v.at(corresponding_point_in_v(u_point_index)); + double m = (projector.x() + projector.y()) / 2.; + return Internal_point(m, m, u_point_index); } inline Internal_point Persistence_graph::get_v_point(int v_point_index) const { - if (!on_the_v_diagonal(v_point_index)) - return v.at(v_point_index); - Internal_point projector = u.at(corresponding_point_in_u(v_point_index)); - double m = (projector.x() + projector.y()) / 2.; - return Internal_point(m,m,v_point_index); + if (!on_the_v_diagonal(v_point_index)) + return v.at(v_point_index); + Internal_point projector = u.at(corresponding_point_in_u(v_point_index)); + double m = (projector.x() + projector.y()) / 2.; + return Internal_point(m, m, v_point_index); } inline double Persistence_graph::diameter_bound() const { - double max = 0.; - for(auto it = u.cbegin(); it != u.cend(); it++) - max = std::max(max, it->y()); - for(auto it = v.cbegin(); it != v.cend(); it++) - max = std::max(max, it->y()); - return max; + double max = 0.; + for (auto it = u.cbegin(); it != u.cend(); it++) + max = std::max(max, it->y()); + for (auto it = v.cbegin(); it != v.cend(); it++) + max = std::max(max, it->y()); + return max; } - } // namespace persistence_diagram } // namespace Gudhi diff --git a/src/Bottleneck_distance/test/CMakeLists.txt b/src/Bottleneck_distance/test/CMakeLists.txt index 13213075..a6979d3c 100644 --- a/src/Bottleneck_distance/test/CMakeLists.txt +++ b/src/Bottleneck_distance/test/CMakeLists.txt @@ -17,7 +17,6 @@ if(CGAL_FOUND) if (NOT CGAL_VERSION VERSION_LESS 4.8.0) if (EIGEN3_FOUND) add_executable ( bottleneckUT bottleneck_unit_test.cpp ) - add_executable ( bottleneck_chrono bottleneck_chrono.cpp ) target_link_libraries(bottleneckUT ${Boost_SYSTEM_LIBRARY} ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY}) # Unitary tests diff --git a/src/Bottleneck_distance/test/bottleneck_chrono.cpp b/src/Bottleneck_distance/test/bottleneck_chrono.cpp deleted file mode 100644 index a30d42b5..00000000 --- a/src/Bottleneck_distance/test/bottleneck_chrono.cpp +++ /dev/null @@ -1,62 +0,0 @@ -/* 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: Francois Godi - * - * Copyright (C) 2015 INRIA (France) - * - * This program is free software: you can redistribute it and/or modify - * it under the terms of the GNU General Public License as published by - * the Free Software Foundation, either version 3 of the License, or - * (at your option) any later version. - * - * This program is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - * GNU General Public License for more details. - * - * You should have received a copy of the GNU General Public License - * along with this program. If not, see . - */ - -#include -#include -#include -#include - -using namespace Gudhi::persistence_diagram; - - -double upper_bound = 400.; // any real >0 - -int main(){ - std::ofstream objetfichier; - objetfichier.open("results.csv", std::ios::out); - - for(int n = 1000; n<=10000; n+=1000){ - std::uniform_real_distribution unif1(0.,upper_bound); - std::uniform_real_distribution unif2(upper_bound/1000.,upper_bound/100.); - std::default_random_engine re; - std::vector< std::pair > v1, v2; - for (int i = 0; i < n; i++) { - double a = unif1(re); - double b = unif1(re); - double x = unif2(re); - double y = unif2(re); - v1.emplace_back(std::min(a,b), std::max(a,b)); - v2.emplace_back(std::min(a,b)+std::min(x,y), std::max(a,b)+std::max(x,y)); - if(i%5==0) - v1.emplace_back(std::min(a,b),std::min(a,b)+x); - if(i%3==0) - v2.emplace_back(std::max(a,b),std::max(a,b)+y); - } - std::chrono::steady_clock::time_point start = std::chrono::steady_clock::now(); - double b = bottleneck_distance(v1, v2); - std::chrono::steady_clock::time_point end = std::chrono::steady_clock::now(); - typedef std::chrono::duration millisecs_t; - millisecs_t duration(std::chrono::duration_cast(end-start)); - objetfichier << n << ";" << duration.count() << ";" << b << std::endl; - } - objetfichier.close(); -} diff --git a/src/Bottleneck_distance/test/bottleneck_unit_test.cpp b/src/Bottleneck_distance/test/bottleneck_unit_test.cpp index fba1d369..e39613b3 100644 --- a/src/Bottleneck_distance/test/bottleneck_unit_test.cpp +++ b/src/Bottleneck_distance/test/bottleneck_unit_test.cpp @@ -4,7 +4,7 @@ * * Author: Francois Godi * - * Copyright (C) 2015 INRIA (France) + * Copyright (C) 2015 INRIA * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by @@ -30,138 +30,138 @@ using namespace Gudhi::persistence_diagram; -int n1 = 81; // a natural number >0 -int n2 = 180; // a natural number >0 -double upper_bound = 406.43; // any real >0 +int n1 = 81; // a natural number >0 +int n2 = 180; // a natural number >0 +double upper_bound = 406.43; // any real >0 -std::uniform_real_distribution unif(0.,upper_bound); +std::uniform_real_distribution unif(0., upper_bound); std::default_random_engine re; std::vector< std::pair > v1, v2; -BOOST_AUTO_TEST_CASE(persistence_graph){ - // Random construction - for (int i = 0; i < n1; i++) { - double a = unif(re); - double b = unif(re); - v1.emplace_back(std::min(a,b), std::max(a,b)); - } - for (int i = 0; i < n2; i++) { - double a = unif(re); - double b = unif(re); - v2.emplace_back(std::min(a,b), std::max(a,b)); - } - Persistence_graph g(v1, v2, 0.); - std::vector d(g.sorted_distances()); - // - BOOST_CHECK(!g.on_the_u_diagonal(n1-1)); - BOOST_CHECK(!g.on_the_u_diagonal(n1)); - BOOST_CHECK(!g.on_the_u_diagonal(n2-1)); - BOOST_CHECK(g.on_the_u_diagonal(n2)); - BOOST_CHECK(!g.on_the_v_diagonal(n1-1)); - BOOST_CHECK(g.on_the_v_diagonal(n1)); - BOOST_CHECK(g.on_the_v_diagonal(n2-1)); - BOOST_CHECK(g.on_the_v_diagonal(n2)); - // - BOOST_CHECK(g.corresponding_point_in_u(0)==n2); - BOOST_CHECK(g.corresponding_point_in_u(n1)==0); - BOOST_CHECK(g.corresponding_point_in_v(0)==n1); - BOOST_CHECK(g.corresponding_point_in_v(n2)==0); - // - BOOST_CHECK(g.size()==(n1+n2)); - // - BOOST_CHECK((int) d.size() == (n1+n2)*(n1+n2) + n1 + n2 + 1); - BOOST_CHECK(std::count(d.begin(), d.end(), g.distance(0,0))>0); - BOOST_CHECK(std::count(d.begin(), d.end(), g.distance(0,n1-1))>0); - BOOST_CHECK(std::count(d.begin(), d.end(), g.distance(0,n1))>0); - BOOST_CHECK(std::count(d.begin(), d.end(), g.distance(0,n2-1))>0); - BOOST_CHECK(std::count(d.begin(), d.end(), g.distance(0,n2))>0); - BOOST_CHECK(std::count(d.begin(), d.end(), g.distance(0,(n1+n2)-1))>0); - BOOST_CHECK(std::count(d.begin(), d.end(), g.distance(n1,0))>0); - BOOST_CHECK(std::count(d.begin(), d.end(), g.distance(n1,n1-1))>0); - BOOST_CHECK(std::count(d.begin(), d.end(), g.distance(n1,n1))>0); - BOOST_CHECK(std::count(d.begin(), d.end(), g.distance(n1,n2-1))>0); - BOOST_CHECK(std::count(d.begin(), d.end(), g.distance(n1,n2))>0); - BOOST_CHECK(std::count(d.begin(), d.end(), g.distance(n1,(n1+n2)-1))>0); - BOOST_CHECK(std::count(d.begin(), d.end(), g.distance((n1+n2)-1,0))>0); - BOOST_CHECK(std::count(d.begin(), d.end(), g.distance((n1+n2)-1,n1-1))>0); - BOOST_CHECK(std::count(d.begin(), d.end(), g.distance((n1+n2)-1,n1))>0); - BOOST_CHECK(std::count(d.begin(), d.end(), g.distance((n1+n2)-1,n2-1))>0); - BOOST_CHECK(std::count(d.begin(), d.end(), g.distance((n1+n2)-1,n2))>0); - BOOST_CHECK(std::count(d.begin(), d.end(), g.distance((n1+n2)-1,(n1+n2)-1))>0); +BOOST_AUTO_TEST_CASE(persistence_graph) { + // Random construction + for (int i = 0; i < n1; i++) { + double a = unif(re); + double b = unif(re); + v1.emplace_back(std::min(a, b), std::max(a, b)); + } + for (int i = 0; i < n2; i++) { + double a = unif(re); + double b = unif(re); + v2.emplace_back(std::min(a, b), std::max(a, b)); + } + Persistence_graph g(v1, v2, 0.); + std::vector d(g.sorted_distances()); + // + BOOST_CHECK(!g.on_the_u_diagonal(n1 - 1)); + BOOST_CHECK(!g.on_the_u_diagonal(n1)); + BOOST_CHECK(!g.on_the_u_diagonal(n2 - 1)); + BOOST_CHECK(g.on_the_u_diagonal(n2)); + BOOST_CHECK(!g.on_the_v_diagonal(n1 - 1)); + BOOST_CHECK(g.on_the_v_diagonal(n1)); + BOOST_CHECK(g.on_the_v_diagonal(n2 - 1)); + BOOST_CHECK(g.on_the_v_diagonal(n2)); + // + BOOST_CHECK(g.corresponding_point_in_u(0) == n2); + BOOST_CHECK(g.corresponding_point_in_u(n1) == 0); + BOOST_CHECK(g.corresponding_point_in_v(0) == n1); + BOOST_CHECK(g.corresponding_point_in_v(n2) == 0); + // + BOOST_CHECK(g.size() == (n1 + n2)); + // + BOOST_CHECK((int) d.size() == (n1 + n2)*(n1 + n2) + n1 + n2 + 1); + BOOST_CHECK(std::count(d.begin(), d.end(), g.distance(0, 0)) > 0); + BOOST_CHECK(std::count(d.begin(), d.end(), g.distance(0, n1 - 1)) > 0); + BOOST_CHECK(std::count(d.begin(), d.end(), g.distance(0, n1)) > 0); + BOOST_CHECK(std::count(d.begin(), d.end(), g.distance(0, n2 - 1)) > 0); + BOOST_CHECK(std::count(d.begin(), d.end(), g.distance(0, n2)) > 0); + BOOST_CHECK(std::count(d.begin(), d.end(), g.distance(0, (n1 + n2) - 1)) > 0); + BOOST_CHECK(std::count(d.begin(), d.end(), g.distance(n1, 0)) > 0); + BOOST_CHECK(std::count(d.begin(), d.end(), g.distance(n1, n1 - 1)) > 0); + BOOST_CHECK(std::count(d.begin(), d.end(), g.distance(n1, n1)) > 0); + BOOST_CHECK(std::count(d.begin(), d.end(), g.distance(n1, n2 - 1)) > 0); + BOOST_CHECK(std::count(d.begin(), d.end(), g.distance(n1, n2)) > 0); + BOOST_CHECK(std::count(d.begin(), d.end(), g.distance(n1, (n1 + n2) - 1)) > 0); + BOOST_CHECK(std::count(d.begin(), d.end(), g.distance((n1 + n2) - 1, 0)) > 0); + BOOST_CHECK(std::count(d.begin(), d.end(), g.distance((n1 + n2) - 1, n1 - 1)) > 0); + BOOST_CHECK(std::count(d.begin(), d.end(), g.distance((n1 + n2) - 1, n1)) > 0); + BOOST_CHECK(std::count(d.begin(), d.end(), g.distance((n1 + n2) - 1, n2 - 1)) > 0); + BOOST_CHECK(std::count(d.begin(), d.end(), g.distance((n1 + n2) - 1, n2)) > 0); + BOOST_CHECK(std::count(d.begin(), d.end(), g.distance((n1 + n2) - 1, (n1 + n2) - 1)) > 0); } BOOST_AUTO_TEST_CASE(neighbors_finder) { - Persistence_graph g(v1, v2, 0.); - Neighbors_finder nf(g, 1.); - for(int v_point_index=1; v_point_index<((n2+n1)*9/10); v_point_index+=2) - nf.add(v_point_index); - // - int v_point_index_1 = nf.pull_near(n2/2); - BOOST_CHECK((v_point_index_1 == -1) || (g.distance(n2/2,v_point_index_1)<=1.)); - std::vector l = nf.pull_all_near(n2/2); - bool v = true; - for(auto it = l.cbegin(); it != l.cend(); ++it) - v = v && (g.distance(n2/2,*it)>1.); - BOOST_CHECK(v); - int v_point_index_2 = nf.pull_near(n2/2); - BOOST_CHECK(v_point_index_2 == -1); + Persistence_graph g(v1, v2, 0.); + Neighbors_finder nf(g, 1.); + for (int v_point_index = 1; v_point_index < ((n2 + n1)*9 / 10); v_point_index += 2) + nf.add(v_point_index); + // + int v_point_index_1 = nf.pull_near(n2 / 2); + BOOST_CHECK((v_point_index_1 == -1) || (g.distance(n2 / 2, v_point_index_1) <= 1.)); + std::vector l = nf.pull_all_near(n2 / 2); + bool v = true; + for (auto it = l.cbegin(); it != l.cend(); ++it) + v = v && (g.distance(n2 / 2, *it) > 1.); + BOOST_CHECK(v); + int v_point_index_2 = nf.pull_near(n2 / 2); + BOOST_CHECK(v_point_index_2 == -1); } BOOST_AUTO_TEST_CASE(layered_neighbors_finder) { - Persistence_graph g(v1, v2, 0.); - Layered_neighbors_finder lnf(g, 1.); - for(int v_point_index=1; v_point_index<((n2+n1)*9/10); v_point_index+=2) - lnf.add(v_point_index, v_point_index % 7); - // - int v_point_index_1 = lnf.pull_near(n2/2,6); - BOOST_CHECK((v_point_index_1 == -1) || (g.distance(n2/2,v_point_index_1)<=1.)); - int v_point_index_2 = lnf.pull_near(n2/2,6); - BOOST_CHECK(v_point_index_2 == -1); - v_point_index_1 = lnf.pull_near(n2/2,0); - BOOST_CHECK((v_point_index_1 == -1) || (g.distance(n2/2,v_point_index_1)<=1.)); - v_point_index_2 = lnf.pull_near(n2/2,0); - BOOST_CHECK(v_point_index_2 == -1); + Persistence_graph g(v1, v2, 0.); + Layered_neighbors_finder lnf(g, 1.); + for (int v_point_index = 1; v_point_index < ((n2 + n1)*9 / 10); v_point_index += 2) + lnf.add(v_point_index, v_point_index % 7); + // + int v_point_index_1 = lnf.pull_near(n2 / 2, 6); + BOOST_CHECK((v_point_index_1 == -1) || (g.distance(n2 / 2, v_point_index_1) <= 1.)); + int v_point_index_2 = lnf.pull_near(n2 / 2, 6); + BOOST_CHECK(v_point_index_2 == -1); + v_point_index_1 = lnf.pull_near(n2 / 2, 0); + BOOST_CHECK((v_point_index_1 == -1) || (g.distance(n2 / 2, v_point_index_1) <= 1.)); + v_point_index_2 = lnf.pull_near(n2 / 2, 0); + BOOST_CHECK(v_point_index_2 == -1); } BOOST_AUTO_TEST_CASE(graph_matching) { - Persistence_graph g(v1, v2, 0.); - Graph_matching m1(g); - m1.set_r(0.); - int e = 0; - while (m1.multi_augment()) - ++e; - BOOST_CHECK(e > 0); - BOOST_CHECK(e <= 2*sqrt(2*(n1+n2))); - Graph_matching m2 = m1; - BOOST_CHECK(!m2.multi_augment()); - m2.set_r(upper_bound); - e = 0; - while (m2.multi_augment()) - ++e; - BOOST_CHECK(e <= 2*sqrt(2*(n1+n2))); - BOOST_CHECK(m2.perfect()); - BOOST_CHECK(!m1.perfect()); + Persistence_graph g(v1, v2, 0.); + Graph_matching m1(g); + m1.set_r(0.); + int e = 0; + while (m1.multi_augment()) + ++e; + BOOST_CHECK(e > 0); + BOOST_CHECK(e <= 2 * sqrt(2 * (n1 + n2))); + Graph_matching m2 = m1; + BOOST_CHECK(!m2.multi_augment()); + m2.set_r(upper_bound); + e = 0; + while (m2.multi_augment()) + ++e; + BOOST_CHECK(e <= 2 * sqrt(2 * (n1 + n2))); + BOOST_CHECK(m2.perfect()); + BOOST_CHECK(!m1.perfect()); } -BOOST_AUTO_TEST_CASE(global){ - std::uniform_real_distribution unif1(0.,upper_bound); - std::uniform_real_distribution unif2(upper_bound/10000.,upper_bound/100.); - std::default_random_engine re; - std::vector< std::pair > v1, v2; - for (int i = 0; i < n1; i++) { - double a = unif1(re); - double b = unif1(re); - double x = unif2(re); - double y = unif2(re); - v1.emplace_back(std::min(a,b), std::max(a,b)); - v2.emplace_back(std::min(a,b)+std::min(x,y), std::max(a,b)+std::max(x,y)); - if(i%5==0) - v1.emplace_back(std::min(a,b),std::min(a,b)+x); - if(i%3==0) - v2.emplace_back(std::max(a,b),std::max(a,b)+y); - } - BOOST_CHECK(bottleneck_distance(v1, v2, 0.) <= upper_bound/100.); - BOOST_CHECK(bottleneck_distance(v1, v2, upper_bound/10000.) <= upper_bound/100. + upper_bound/10000.); - BOOST_CHECK(std::abs(bottleneck_distance(v1, v2, 0.) - bottleneck_distance(v1, v2, upper_bound/10000.)) <= upper_bound/10000.); +BOOST_AUTO_TEST_CASE(global) { + std::uniform_real_distribution unif1(0., upper_bound); + std::uniform_real_distribution unif2(upper_bound / 10000., upper_bound / 100.); + std::default_random_engine re; + std::vector< std::pair > v1, v2; + for (int i = 0; i < n1; i++) { + double a = unif1(re); + double b = unif1(re); + double x = unif2(re); + double y = unif2(re); + v1.emplace_back(std::min(a, b), std::max(a, b)); + v2.emplace_back(std::min(a, b) + std::min(x, y), std::max(a, b) + std::max(x, y)); + if (i % 5 == 0) + v1.emplace_back(std::min(a, b), std::min(a, b) + x); + if (i % 3 == 0) + v2.emplace_back(std::max(a, b), std::max(a, b) + y); + } + BOOST_CHECK(bottleneck_distance(v1, v2, 0.) <= upper_bound / 100.); + BOOST_CHECK(bottleneck_distance(v1, v2, upper_bound / 10000.) <= upper_bound / 100. + upper_bound / 10000.); + BOOST_CHECK(std::abs(bottleneck_distance(v1, v2, 0.) - bottleneck_distance(v1, v2, upper_bound / 10000.)) <= upper_bound / 10000.); } diff --git a/src/Tangential_complex/benchmark/CMakeLists.txt b/src/Tangential_complex/benchmark/CMakeLists.txt index 56dd8128..788c2b4d 100644 --- a/src/Tangential_complex/benchmark/CMakeLists.txt +++ b/src/Tangential_complex/benchmark/CMakeLists.txt @@ -1,15 +1,6 @@ cmake_minimum_required(VERSION 2.6) project(Tangential_complex_benchmark) -if (GCOVR_PATH) - # for gcovr to make coverage reports - Corbera Jenkins plugin - set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -fprofile-arcs -ftest-coverage") -endif() -if (GPROF_PATH) - # for gprof to make coverage reports - Jenkins - set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -pg") -endif() - # need CGAL 4.8 if(CGAL_FOUND) if (NOT CGAL_VERSION VERSION_LESS 4.8.0) diff --git a/src/common/include/gudhi_patches/Bottleneck_distance_CGAL_patches.txt b/src/common/include/gudhi_patches/Bottleneck_distance_CGAL_patches.txt new file mode 100644 index 00000000..a588d113 --- /dev/null +++ b/src/common/include/gudhi_patches/Bottleneck_distance_CGAL_patches.txt @@ -0,0 +1,3 @@ +CGAL/Kd_tree.h +CGAL/Kd_tree_node.h +CGAL/Orthogonal_incremental_neighbor_search.h diff --git a/src/common/include/gudhi_patches/CGAL/Kd_tree.h b/src/common/include/gudhi_patches/CGAL/Kd_tree.h new file mode 100644 index 00000000..f085b0da --- /dev/null +++ b/src/common/include/gudhi_patches/CGAL/Kd_tree.h @@ -0,0 +1,582 @@ +// Copyright (c) 2002,2011,2014 Utrecht University (The Netherlands), Max-Planck-Institute Saarbruecken (Germany). +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org). +// You can redistribute it and/or modify it under the terms of the GNU +// General Public License as published by the Free Software Foundation, +// either version 3 of the License, or (at your option) any later version. +// +// Licensees holding a valid commercial license may use this file in +// accordance with the commercial license agreement provided with the software. +// +// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE +// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +// +// $URL$ +// $Id$ +// +// Author(s) : Hans Tangelder (), +// : Waqar Khan + +#ifndef CGAL_KD_TREE_H +#define CGAL_KD_TREE_H + +#include "Kd_tree_node.h" + +#include +#include +#include + +#include +#include +#include + + +#include +#include +#include + +#ifdef CGAL_HAS_THREADS +#include +#endif + +namespace CGAL { + +//template , class UseExtendedNode = Tag_true > +template , class UseExtendedNode = Tag_true > +class Kd_tree { + +public: + typedef SearchTraits Traits; + typedef Splitter_ Splitter; + typedef typename SearchTraits::Point_d Point_d; + typedef typename Splitter::Container Point_container; + + typedef typename SearchTraits::FT FT; + typedef Kd_tree_node Node; + typedef Kd_tree_leaf_node Leaf_node; + typedef Kd_tree_internal_node Internal_node; + typedef Kd_tree Tree; + typedef Kd_tree Self; + + typedef Node* Node_handle; + typedef const Node* Node_const_handle; + typedef Leaf_node* Leaf_node_handle; + typedef const Leaf_node* Leaf_node_const_handle; + typedef Internal_node* Internal_node_handle; + typedef const Internal_node* Internal_node_const_handle; + typedef typename std::vector::const_iterator Point_d_iterator; + typedef typename std::vector::const_iterator Point_d_const_iterator; + typedef typename Splitter::Separator Separator; + typedef typename std::vector::const_iterator iterator; + typedef typename std::vector::const_iterator const_iterator; + + typedef typename std::vector::size_type size_type; + + typedef typename internal::Get_dimension_tag::Dimension D; + +private: + SearchTraits traits_; + Splitter split; + + + // wokaround for https://svn.boost.org/trac/boost/ticket/9332 +#if (_MSC_VER == 1800) && (BOOST_VERSION == 105500) + std::deque internal_nodes; + std::deque leaf_nodes; +#else + boost::container::deque internal_nodes; + boost::container::deque leaf_nodes; +#endif + + Node_handle tree_root; + + Kd_tree_rectangle* bbox; + std::vector pts; + + // Instead of storing the points in arrays in the Kd_tree_node + // we put all the data in a vector in the Kd_tree. + // and we only store an iterator range in the Kd_tree_node. + // + std::vector data; + + + #ifdef CGAL_HAS_THREADS + mutable CGAL_MUTEX building_mutex;//mutex used to protect const calls inducing build() + #endif + bool built_; + bool removed_; + + // protected copy constructor + Kd_tree(const Tree& tree) + : traits_(tree.traits_),built_(tree.built_) + {}; + + + // Instead of the recursive construction of the tree in the class Kd_tree_node + // we do this in the tree class. The advantage is that we then can optimize + // the allocation of the nodes. + + // The leaf node + Node_handle + create_leaf_node(Point_container& c) + { + Leaf_node node(true , static_cast(c.size())); + std::ptrdiff_t tmp = c.begin() - data.begin(); + node.data = pts.begin() + tmp; + + leaf_nodes.push_back(node); + Leaf_node_handle nh = &leaf_nodes.back(); + + + return nh; + } + + + // The internal node + + Node_handle + create_internal_node(Point_container& c, const Tag_true&) + { + return create_internal_node_use_extension(c); + } + + Node_handle + create_internal_node(Point_container& c, const Tag_false&) + { + return create_internal_node(c); + } + + + + // TODO: Similiar to the leaf_init function above, a part of the code should be + // moved to a the class Kd_tree_node. + // It is not proper yet, but the goal was to see if there is + // a potential performance gain through the Compact_container + Node_handle + create_internal_node_use_extension(Point_container& c) + { + Internal_node node(false); + internal_nodes.push_back(node); + Internal_node_handle nh = &internal_nodes.back(); + + Separator sep; + Point_container c_low(c.dimension(),traits_); + split(sep, c, c_low); + nh->set_separator(sep); + + int cd = nh->cutting_dimension(); + if(!c_low.empty()){ + nh->lower_low_val = c_low.tight_bounding_box().min_coord(cd); + nh->lower_high_val = c_low.tight_bounding_box().max_coord(cd); + } + else{ + nh->lower_low_val = nh->cutting_value(); + nh->lower_high_val = nh->cutting_value(); + } + if(!c.empty()){ + nh->upper_low_val = c.tight_bounding_box().min_coord(cd); + nh->upper_high_val = c.tight_bounding_box().max_coord(cd); + } + else{ + nh->upper_low_val = nh->cutting_value(); + nh->upper_high_val = nh->cutting_value(); + } + + CGAL_assertion(nh->cutting_value() >= nh->lower_low_val); + CGAL_assertion(nh->cutting_value() <= nh->upper_high_val); + + if (c_low.size() > split.bucket_size()){ + nh->lower_ch = create_internal_node_use_extension(c_low); + }else{ + nh->lower_ch = create_leaf_node(c_low); + } + if (c.size() > split.bucket_size()){ + nh->upper_ch = create_internal_node_use_extension(c); + }else{ + nh->upper_ch = create_leaf_node(c); + } + + + + + return nh; + } + + + // Note also that I duplicated the code to get rid if the if's for + // the boolean use_extension which was constant over the construction + Node_handle + create_internal_node(Point_container& c) + { + Internal_node node(false); + internal_nodes.push_back(node); + Internal_node_handle nh = &internal_nodes.back(); + Separator sep; + + Point_container c_low(c.dimension(),traits_); + split(sep, c, c_low); + nh->set_separator(sep); + + if (c_low.size() > split.bucket_size()){ + nh->lower_ch = create_internal_node(c_low); + }else{ + nh->lower_ch = create_leaf_node(c_low); + } + if (c.size() > split.bucket_size()){ + nh->upper_ch = create_internal_node(c); + }else{ + nh->upper_ch = create_leaf_node(c); + } + + + + return nh; + } + + + +public: + + Kd_tree(Splitter s = Splitter(),const SearchTraits traits=SearchTraits()) + : traits_(traits),split(s), built_(false), removed_(false) + {} + + template + Kd_tree(InputIterator first, InputIterator beyond, + Splitter s = Splitter(),const SearchTraits traits=SearchTraits()) + : traits_(traits),split(s), built_(false), removed_(false) + { + pts.insert(pts.end(), first, beyond); + } + + bool empty() const { + return pts.empty(); + } + + void + build() + { + // This function is not ready to be called when a tree already exists, one + // must call invalidate_built() first. + CGAL_assertion(!is_built()); + CGAL_assertion(!removed_); + const Point_d& p = *pts.begin(); + typename SearchTraits::Construct_cartesian_const_iterator_d ccci=traits_.construct_cartesian_const_iterator_d_object(); + int dim = static_cast(std::distance(ccci(p), ccci(p,0))); + + data.reserve(pts.size()); + for(unsigned int i = 0; i < pts.size(); i++){ + data.push_back(&pts[i]); + } + Point_container c(dim, data.begin(), data.end(),traits_); + bbox = new Kd_tree_rectangle(c.bounding_box()); + if (c.size() <= split.bucket_size()){ + tree_root = create_leaf_node(c); + }else { + tree_root = create_internal_node(c, UseExtendedNode()); + } + + //Reorder vector for spatial locality + std::vector ptstmp; + ptstmp.resize(pts.size()); + for (std::size_t i = 0; i < pts.size(); ++i){ + ptstmp[i] = *data[i]; + } + for(std::size_t i = 0; i < leaf_nodes.size(); ++i){ + std::ptrdiff_t tmp = leaf_nodes[i].begin() - pts.begin(); + leaf_nodes[i].data = ptstmp.begin() + tmp; + } + pts.swap(ptstmp); + + data.clear(); + + built_ = true; + } + +private: + //any call to this function is for the moment not threadsafe + void const_build() const { + #ifdef CGAL_HAS_THREADS + //this ensure that build() will be called once + CGAL_SCOPED_LOCK(building_mutex); + if(!is_built()) + #endif + const_cast(this)->build(); //THIS IS NOT THREADSAFE + } +public: + + bool is_built() const + { + return built_; + } + + void invalidate_built() + { + if(removed_){ + // Walk the tree to collect the remaining points. + // Writing directly to pts would likely work, but better be safe. + std::vector ptstmp; + //ptstmp.resize(root()->num_items()); + root()->tree_items(std::back_inserter(ptstmp)); + pts.swap(ptstmp); + removed_=false; + CGAL_assertion(is_built()); // the rest of the cleanup must happen + } + if(is_built()){ + internal_nodes.clear(); + leaf_nodes.clear(); + data.clear(); + delete bbox; + built_ = false; + } + } + + void clear() + { + invalidate_built(); + pts.clear(); + removed_ = false; + } + + void + insert(const Point_d& p) + { + invalidate_built(); + pts.push_back(p); + } + + template + void + insert(InputIterator first, InputIterator beyond) + { + invalidate_built(); + pts.insert(pts.end(),first, beyond); + } + +private: + struct Equal_by_coordinates { + SearchTraits const* traits; + Point_d const* pp; + bool operator()(Point_d const&q) const { + typename SearchTraits::Construct_cartesian_const_iterator_d ccci=traits->construct_cartesian_const_iterator_d_object(); + return std::equal(ccci(*pp), ccci(*pp,0), ccci(q)); + } + }; + Equal_by_coordinates equal_by_coordinates(Point_d const&p){ + Equal_by_coordinates ret = { &traits(), &p }; + return ret; + } + +public: + void + remove(const Point_d& p) + { + remove(p, equal_by_coordinates(p)); + } + + template + void + remove(const Point_d& p, Equal const& equal_to_p) + { +#if 0 + // This code could have quadratic runtime. + if (!is_built()) { + std::vector::iterator pi = std::find(pts.begin(), pts.end(), p); + // Precondition: the point must be there. + CGAL_assertion (pi != pts.end()); + pts.erase(pi); + return; + } +#endif + bool success = remove_(p, 0, false, 0, false, root(), equal_to_p); + CGAL_assertion(success); + + // Do not set the flag is the tree has been cleared. + if(is_built()) + removed_ |= success; + } +private: + template + bool remove_(const Point_d& p, + Internal_node_handle grandparent, bool parent_islower, + Internal_node_handle parent, bool islower, + Node_handle node, Equal const& equal_to_p) { + // Recurse to locate the point + if (!node->is_leaf()) { + Internal_node_handle newparent = static_cast(node); + // FIXME: This should be if(xcutting_dimension()] <= newparent->cutting_value()) { + if (remove_(p, parent, islower, newparent, true, newparent->lower(), equal_to_p)) + return true; + } + //if (traits().construct_cartesian_const_iterator_d_object()(p)[newparent->cutting_dimension()] >= newparent->cutting_value()) + return remove_(p, parent, islower, newparent, false, newparent->upper(), equal_to_p); + + CGAL_assertion(false); // Point was not found + } + + // Actual removal + Leaf_node_handle lnode = static_cast(node); + if (lnode->size() > 1) { + iterator pi = std::find_if(lnode->begin(), lnode->end(), equal_to_p); + // FIXME: we should ensure this never happens + if (pi == lnode->end()) return false; + iterator lasti = lnode->end() - 1; + if (pi != lasti) { + // Hack to get a non-const iterator + std::iter_swap(pts.begin()+(pi-pts.begin()), pts.begin()+(lasti-pts.begin())); + } + lnode->drop_last_point(); + } else if (!equal_to_p(*lnode->begin())) { + // FIXME: we should ensure this never happens + return false; + } else if (grandparent) { + Node_handle brother = islower ? parent->upper() : parent->lower(); + if (parent_islower) + grandparent->set_lower(brother); + else + grandparent->set_upper(brother); + } else if (parent) { + tree_root = islower ? parent->upper() : parent->lower(); + } else { + clear(); + } + return true; + } + +public: + //For efficiency; reserve the size of the points vectors in advance (if the number of points is already known). + void reserve(size_t size) + { + pts.reserve(size); + } + + //Get the capacity of the underlying points vector. + size_t capacity() + { + return pts.capacity(); + } + + + template + OutputIterator + search(OutputIterator it, const FuzzyQueryItem& q) const + { + if(! pts.empty()){ + + if(! is_built()){ + const_build(); + } + Kd_tree_rectangle b(*bbox); + return tree_root->search(it,q,b); + } + return it; + } + + + template + boost::optional + search_any_point(const FuzzyQueryItem& q) const + { + if(! pts.empty()){ + + if(! is_built()){ + const_build(); + } + Kd_tree_rectangle b(*bbox); + return tree_root->search_any_point(q,b); + } + return boost::none; + } + + + ~Kd_tree() { + if(is_built()){ + delete bbox; + } + } + + + const SearchTraits& + traits() const + { + return traits_; + } + + Node_const_handle + root() const + { + if(! is_built()){ + const_build(); + } + return tree_root; + } + + Node_handle + root() + { + if(! is_built()){ + build(); + } + return tree_root; + } + + void + print() const + { + if(! is_built()){ + const_build(); + } + root()->print(); + } + + const Kd_tree_rectangle& + bounding_box() const + { + if(! is_built()){ + const_build(); + } + return *bbox; + } + + const_iterator + begin() const + { + return pts.begin(); + } + + const_iterator + end() const + { + return pts.end(); + } + + size_type + size() const + { + return pts.size(); + } + + // Print statistics of the tree. + std::ostream& + statistics(std::ostream& s) const + { + if(! is_built()){ + const_build(); + } + s << "Tree statistics:" << std::endl; + s << "Number of items stored: " + << root()->num_items() << std::endl; + s << "Number of nodes: " + << root()->num_nodes() << std::endl; + s << " Tree depth: " << root()->depth() << std::endl; + return s; + } + + +}; + +} // namespace CGAL + +#endif // CGAL_KD_TREE_H diff --git a/src/common/include/gudhi_patches/CGAL/Kd_tree_node.h b/src/common/include/gudhi_patches/CGAL/Kd_tree_node.h new file mode 100644 index 00000000..909ee260 --- /dev/null +++ b/src/common/include/gudhi_patches/CGAL/Kd_tree_node.h @@ -0,0 +1,586 @@ +// Copyright (c) 2002,2011 Utrecht University (The Netherlands). +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org). +// You can redistribute it and/or modify it under the terms of the GNU +// General Public License as published by the Free Software Foundation, +// either version 3 of the License, or (at your option) any later version. +// +// Licensees holding a valid commercial license may use this file in +// accordance with the commercial license agreement provided with the software. +// +// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE +// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +// +// $URL$ +// $Id$ +// +// +// Authors : Hans Tangelder () + +#ifndef CGAL_KD_TREE_NODE_H +#define CGAL_KD_TREE_NODE_H + +#include "CGAL/Splitters.h" + +#include +#include + +namespace CGAL { + + template + class Kd_tree; + + template < class TreeTraits, class Splitter, class UseExtendedNode > + class Kd_tree_node { + + friend class Kd_tree; + + typedef typename Kd_tree::Node_handle Node_handle; + typedef typename Kd_tree::Node_const_handle Node_const_handle; + typedef typename Kd_tree::Internal_node_handle Internal_node_handle; + typedef typename Kd_tree::Internal_node_const_handle Internal_node_const_handle; + typedef typename Kd_tree::Leaf_node_handle Leaf_node_handle; + typedef typename Kd_tree::Leaf_node_const_handle Leaf_node_const_handle; + typedef typename TreeTraits::Point_d Point_d; + + typedef typename TreeTraits::FT FT; + typedef typename Kd_tree::Separator Separator; + typedef typename Kd_tree::Point_d_iterator Point_d_iterator; + typedef typename Kd_tree::iterator iterator; + typedef typename Kd_tree::D D; + + bool leaf; + + public : + Kd_tree_node(bool leaf_) + :leaf(leaf_){} + + bool is_leaf() const{ + return leaf; + } + + std::size_t + num_items() const + { + if (is_leaf()){ + Leaf_node_const_handle node = + static_cast(this); + return node->size(); + } + else { + Internal_node_const_handle node = + static_cast(this); + return node->lower()->num_items() + node->upper()->num_items(); + } + } + + std::size_t + num_nodes() const + { + if (is_leaf()) return 1; + else { + Internal_node_const_handle node = + static_cast(this); + return node->lower()->num_nodes() + node->upper()->num_nodes(); + } + } + + int + depth(const int current_max_depth) const + { + if (is_leaf()){ + return current_max_depth; + } + else { + Internal_node_const_handle node = + static_cast(this); + return + (std::max)( node->lower()->depth(current_max_depth + 1), + node->upper()->depth(current_max_depth + 1)); + } + } + + int + depth() const + { + return depth(1); + } + + template + OutputIterator + tree_items(OutputIterator it) const { + if (is_leaf()) { + Leaf_node_const_handle node = + static_cast(this); + if (node->size()>0) + for (iterator i=node->begin(); i != node->end(); i++) + {*it=*i; ++it;} + } + else { + Internal_node_const_handle node = + static_cast(this); + it=node->lower()->tree_items(it); + it=node->upper()->tree_items(it); + } + return it; + } + + + boost::optional + any_tree_item() const { + boost::optional result = boost::none; + if (is_leaf()) { + Leaf_node_const_handle node = + static_cast(this); + if (node->size()>0){ + return boost::make_optional(*(node->begin())); + } + } + else { + Internal_node_const_handle node = + static_cast(this); + result = node->lower()->any_tree_item(); + if(! result){ + result = node->upper()->any_tree_item(); + } + } + return result; + } + + + void + indent(int d) const + { + for(int i = 0; i < d; i++){ + std::cout << " "; + } + } + + + void + print(int d = 0) const + { + if (is_leaf()) { + Leaf_node_const_handle node = + static_cast(this); + indent(d); + std::cout << "leaf" << std::endl; + if (node->size()>0) + for (iterator i=node->begin(); i != node->end(); i++) + {indent(d);std::cout << *i << std::endl;} + } + else { + Internal_node_const_handle node = + static_cast(this); + indent(d); + std::cout << "lower tree" << std::endl; + node->lower()->print(d+1); + indent(d); + std::cout << "separator: dim = " << node->cutting_dimension() << " val = " << node->cutting_value() << std::endl; + indent(d); + std::cout << "upper tree" << std::endl; + node->upper()->print(d+1); + } + } + + + template + OutputIterator + search(OutputIterator it, const FuzzyQueryItem& q, + Kd_tree_rectangle& b) const + { + if (is_leaf()) { + Leaf_node_const_handle node = + static_cast(this); + if (node->size()>0) + for (iterator i=node->begin(); i != node->end(); i++) + if (q.contains(*i)) + {*it++=*i;} + } + else { + Internal_node_const_handle node = + static_cast(this); + // after splitting b denotes the lower part of b + Kd_tree_rectangle b_upper(b); + b.split(b_upper, node->cutting_dimension(), + node->cutting_value()); + + if (q.outer_range_contains(b)) + it=node->lower()->tree_items(it); + else + if (q.inner_range_intersects(b)) + it=node->lower()->search(it,q,b); + if (q.outer_range_contains(b_upper)) + it=node->upper()->tree_items(it); + else + if (q.inner_range_intersects(b_upper)) + it=node->upper()->search(it,q,b_upper); + }; + return it; + } + + + template + boost::optional + search_any_point(const FuzzyQueryItem& q, + Kd_tree_rectangle& b) const + { + boost::optional result = boost::none; + if (is_leaf()) { + Leaf_node_const_handle node = + static_cast(this); + if (node->size()>0) + for (iterator i=node->begin(); i != node->end(); i++) + if (q.contains(*i)) + { result = *i; break; } + } + else { + Internal_node_const_handle node = + static_cast(this); + // after splitting b denotes the lower part of b + Kd_tree_rectangle b_upper(b); + b.split(b_upper, node->cutting_dimension(), + node->cutting_value()); + + if (q.outer_range_contains(b)){ + result = node->lower()->any_tree_item(); + }else{ + if (q.inner_range_intersects(b)){ + result = node->lower()->search_any_point(q,b); + } + } + if(result){ + return result; + } + if (q.outer_range_contains(b_upper)){ + result = node->upper()->any_tree_item(); + }else{ + if (q.inner_range_intersects(b_upper)) + result = node->upper()->search_any_point(q,b_upper); + } + } + return result; + } + + }; + + + template < class TreeTraits, class Splitter, class UseExtendedNode > + class Kd_tree_leaf_node : public Kd_tree_node< TreeTraits, Splitter, UseExtendedNode >{ + + friend class Kd_tree; + + typedef typename Kd_tree::iterator iterator; + typedef Kd_tree_node< TreeTraits, Splitter, UseExtendedNode> Base; + typedef typename TreeTraits::Point_d Point_d; + + private: + + // private variables for leaf nodes + boost::int32_t n; // denotes number of items in a leaf node + iterator data; // iterator to data in leaf node + + + public: + + // default constructor + Kd_tree_leaf_node() + {} + + Kd_tree_leaf_node(bool leaf_ ) + : Base(leaf_) + {} + + Kd_tree_leaf_node(bool leaf_,unsigned int n_ ) + : Base(leaf_), n(n_) + {} + + // members for all nodes + + // members for leaf nodes only + inline + unsigned int + size() const + { + return n; + } + + inline + iterator + begin() const + { + return data; + } + + inline + iterator + end() const + { + return data + n; + } + + inline + void + drop_last_point() + { + --n; + } + + }; //leaf node + + + + template < class TreeTraits, class Splitter, class UseExtendedNode> + class Kd_tree_internal_node : public Kd_tree_node< TreeTraits, Splitter, UseExtendedNode >{ + + friend class Kd_tree; + + typedef Kd_tree_node< TreeTraits, Splitter, UseExtendedNode> Base; + typedef typename Kd_tree::Node_handle Node_handle; + typedef typename Kd_tree::Node_const_handle Node_const_handle; + + typedef typename TreeTraits::FT FT; + typedef typename Kd_tree::Separator Separator; + + private: + + // private variables for internal nodes + boost::int32_t cut_dim; + FT cut_val; + Node_handle lower_ch, upper_ch; + + + // private variables for extended internal nodes + FT upper_low_val; + FT upper_high_val; + FT lower_low_val; + FT lower_high_val; + + + public: + + // default constructor + Kd_tree_internal_node() + {} + + Kd_tree_internal_node(bool leaf_) + : Base(leaf_) + {} + + + // members for internal node and extended internal node + + inline + Node_const_handle + lower() const + { + return lower_ch; + } + + inline + Node_const_handle + upper() const + { + return upper_ch; + } + + inline + Node_handle + lower() + { + return lower_ch; + } + + inline + Node_handle + upper() + { + return upper_ch; + } + + inline + void + set_lower(Node_handle nh) + { + lower_ch = nh; + } + + inline + void + set_upper(Node_handle nh) + { + upper_ch = nh; + } + + // inline Separator& separator() {return sep; } + // use instead + inline + void set_separator(Separator& sep){ + cut_dim = sep.cutting_dimension(); + cut_val = sep.cutting_value(); + } + + inline + FT + cutting_value() const + { + return cut_val; + } + + inline + int + cutting_dimension() const + { + return cut_dim; + } + + // members for extended internal node only + inline + FT + upper_low_value() const + { + return upper_low_val; + } + + inline + FT + upper_high_value() const + { + return upper_high_val; + } + + inline + FT + lower_low_value() const + { + return lower_low_val; + } + + inline + FT + lower_high_value() const + { + return lower_high_val; + } + + /*Separator& + separator() + { + return Separator(cutting_dimension,cutting_value); + }*/ + + + };//internal node + + template < class TreeTraits, class Splitter> + class Kd_tree_internal_node : public Kd_tree_node< TreeTraits, Splitter, Tag_false >{ + + friend class Kd_tree; + + typedef Kd_tree_node< TreeTraits, Splitter, Tag_false> Base; + typedef typename Kd_tree::Node_handle Node_handle; + typedef typename Kd_tree::Node_const_handle Node_const_handle; + + typedef typename TreeTraits::FT FT; + typedef typename Kd_tree::Separator Separator; + + private: + + // private variables for internal nodes + boost::uint8_t cut_dim; + FT cut_val; + + Node_handle lower_ch, upper_ch; + + public: + + // default constructor + Kd_tree_internal_node() + {} + + Kd_tree_internal_node(bool leaf_) + : Base(leaf_) + {} + + + // members for internal node and extended internal node + + inline + Node_const_handle + lower() const + { + return lower_ch; + } + + inline + Node_const_handle + upper() const + { + return upper_ch; + } + + inline + Node_handle + lower() + { + return lower_ch; + } + + inline + Node_handle + upper() + { + return upper_ch; + } + + inline + void + set_lower(Node_handle nh) + { + lower_ch = nh; + } + + inline + void + set_upper(Node_handle nh) + { + upper_ch = nh; + } + + // inline Separator& separator() {return sep; } + // use instead + + inline + void set_separator(Separator& sep){ + cut_dim = sep.cutting_dimension(); + cut_val = sep.cutting_value(); + } + + inline + FT + cutting_value() const + { + return cut_val; + } + + inline + int + cutting_dimension() const + { + return cut_dim; + } + + /* Separator& + separator() + { + return Separator(cutting_dimension,cutting_value); + }*/ + + + };//internal node + + + +} // namespace CGAL +#endif // CGAL_KDTREE_NODE_H diff --git a/src/common/include/gudhi_patches/CGAL/Orthogonal_incremental_neighbor_search.h b/src/common/include/gudhi_patches/CGAL/Orthogonal_incremental_neighbor_search.h new file mode 100644 index 00000000..dbe707ed --- /dev/null +++ b/src/common/include/gudhi_patches/CGAL/Orthogonal_incremental_neighbor_search.h @@ -0,0 +1,621 @@ +// Copyright (c) 2002,2011 Utrecht University (The Netherlands). +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org). +// You can redistribute it and/or modify it under the terms of the GNU +// General Public License as published by the Free Software Foundation, +// either version 3 of the License, or (at your option) any later version. +// +// Licensees holding a valid commercial license may use this file in +// accordance with the commercial license agreement provided with the software. +// +// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE +// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +// +// $URL$ +// $Id$ +// +// +// Author(s) : Hans Tangelder () + +#ifndef CGAL_ORTHOGONAL_INCREMENTAL_NEIGHBOR_SEARCH +#define CGAL_ORTHOGONAL_INCREMENTAL_NEIGHBOR_SEARCH + +#include "CGAL/Kd_tree.h" + +#include +#include +#include +#include +#include +#include + +namespace CGAL { + + template ::type, + class Splitter_ = Sliding_midpoint, + class Tree_= Kd_tree > + class Orthogonal_incremental_neighbor_search { + + public: + typedef Splitter_ Splitter; + typedef Tree_ Tree; + typedef Distance_ Distance; + typedef typename SearchTraits::Point_d Point_d; + typedef typename Distance::Query_item Query_item; + typedef typename SearchTraits::FT FT; + typedef typename Tree::Point_d_iterator Point_d_iterator; + typedef typename Tree::Node_const_handle Node_const_handle; + + typedef std::pair Point_with_transformed_distance; + typedef CGAL::cpp11::tuple > Node_with_distance; + typedef std::vector Node_with_distance_vector; + typedef std::vector Point_with_transformed_distance_vector; + + template + struct Object_wrapper + { + T object; + Object_wrapper(const T& t):object(t){} + const T& operator* () const { return object; } + const T* operator-> () const { return &object; } + }; + + class Iterator_implementation { + SearchTraits traits; + public: + + int number_of_neighbours_computed; + int number_of_internal_nodes_visited; + int number_of_leaf_nodes_visited; + int number_of_items_visited; + + private: + + typedef std::vector Distance_vector; + + Distance_vector dists; + + Distance Orthogonal_distance_instance; + + FT multiplication_factor; + + Query_item query_point; + + FT distance_to_root; + + bool search_nearest_neighbour; + + FT rd; + + + class Priority_higher { + public: + + bool search_nearest; + + Priority_higher(bool search_the_nearest_neighbour) + : search_nearest(search_the_nearest_neighbour) + {} + + //highest priority is smallest distance + bool + operator() (Node_with_distance* n1, Node_with_distance* n2) const + { + return (search_nearest) ? (CGAL::cpp11::get<1>(*n1) > CGAL::cpp11::get<1>(*n2)) : (CGAL::cpp11::get<1>(*n2) > CGAL::cpp11::get<1>(*n1)); + } + }; + + class Distance_smaller { + + public: + + bool search_nearest; + + Distance_smaller(bool search_the_nearest_neighbour) + : search_nearest(search_the_nearest_neighbour) + {} + + //highest priority is smallest distance + bool operator() (Point_with_transformed_distance* p1, Point_with_transformed_distance* p2) const + { + return (search_nearest) ? (p1->second > p2->second) : (p2->second > p1->second); + } + }; + + + std::priority_queue PriorityQueue; + + public: + std::priority_queue Item_PriorityQueue; + + + public: + + int reference_count; + + + + // constructor + Iterator_implementation(const Tree& tree,const Query_item& q, const Distance& tr, + FT Eps=FT(0.0), bool search_nearest=true) + : traits(tree.traits()),number_of_neighbours_computed(0), number_of_internal_nodes_visited(0), + number_of_leaf_nodes_visited(0), number_of_items_visited(0), + Orthogonal_distance_instance(tr), multiplication_factor(Orthogonal_distance_instance.transformed_distance(FT(1.0)+Eps)), + query_point(q), search_nearest_neighbour(search_nearest), + PriorityQueue(Priority_higher(search_nearest)), Item_PriorityQueue(Distance_smaller(search_nearest)), + reference_count(1) + + + { + if (tree.empty()) return; + + typename SearchTraits::Construct_cartesian_const_iterator_d ccci=traits.construct_cartesian_const_iterator_d_object(); + int dim = static_cast(std::distance(ccci(q), ccci(q,0))); + + dists.resize(dim); + for(int i=0 ; i(*The_Root); + Compute_the_next_nearest_neighbour(); + } + else{ + distance_to_root= + Orthogonal_distance_instance.max_distance_to_rectangle(q, + tree.bounding_box(), dists); + Node_with_distance *The_Root = new Node_with_distance(tree.root(), + distance_to_root, dists); + PriorityQueue.push(The_Root); + + // rd is the distance of the top of the priority queue to q + rd=CGAL::cpp11::get<1>(*The_Root); + Compute_the_next_furthest_neighbour(); + } + + + } + + // * operator + const Point_with_transformed_distance& + operator* () const + { + return *(Item_PriorityQueue.top()); + } + + // prefix operator + Iterator_implementation& + operator++() + { + Delete_the_current_item_top(); + if(search_nearest_neighbour) + Compute_the_next_nearest_neighbour(); + else + Compute_the_next_furthest_neighbour(); + return *this; + } + + // postfix operator + Object_wrapper + operator++(int) + { + Object_wrapper result( *(Item_PriorityQueue.top()) ); + ++*this; + return result; + } + + // Print statistics of the general priority search process. + std::ostream& + statistics (std::ostream& s) const { + s << "Orthogonal priority search statistics:" + << std::endl; + s << "Number of internal nodes visited:" + << number_of_internal_nodes_visited << std::endl; + s << "Number of leaf nodes visited:" + << number_of_leaf_nodes_visited << std::endl; + s << "Number of items visited:" + << number_of_items_visited << std::endl; + s << "Number of neighbours computed:" + << number_of_neighbours_computed << std::endl; + return s; + } + + + //destructor + ~Iterator_implementation() + { + while (!PriorityQueue.empty()) { + Node_with_distance* The_top=PriorityQueue.top(); + PriorityQueue.pop(); + delete The_top; + } + while (!Item_PriorityQueue.empty()) { + Point_with_transformed_distance* The_top=Item_PriorityQueue.top(); + Item_PriorityQueue.pop(); + delete The_top; + } + } + + private: + + void + Delete_the_current_item_top() + { + Point_with_transformed_distance* The_item_top=Item_PriorityQueue.top(); + Item_PriorityQueue.pop(); + delete The_item_top; + } + + void + Compute_the_next_nearest_neighbour() + { + // compute the next item + bool next_neighbour_found=false; + if (!(Item_PriorityQueue.empty())) { + next_neighbour_found= + (multiplication_factor*rd > Item_PriorityQueue.top()->second); + } + typename SearchTraits::Construct_cartesian_const_iterator_d construct_it=traits.construct_cartesian_const_iterator_d_object(); + typename SearchTraits::Cartesian_const_iterator_d query_point_it = construct_it(query_point); + // otherwise browse the tree further + while ((!next_neighbour_found) && (!PriorityQueue.empty())) { + Node_with_distance* The_node_top=PriorityQueue.top(); + Node_const_handle N= CGAL::cpp11::get<0>(*The_node_top); + dists = CGAL::cpp11::get<2>(*The_node_top); + PriorityQueue.pop(); + delete The_node_top; + FT copy_rd=rd; + while (!(N->is_leaf())) { // compute new distance + typename Tree::Internal_node_const_handle node = + static_cast(N); + number_of_internal_nodes_visited++; + int new_cut_dim=node->cutting_dimension(); + FT new_rd,dst = dists[new_cut_dim]; + FT val = *(query_point_it + new_cut_dim); + FT diff1 = val - node->upper_low_value(); + FT diff2 = val - node->lower_high_value(); + if (diff1 + diff2 < FT(0.0)) { + new_rd= + Orthogonal_distance_instance.new_distance(copy_rd,dst,diff1,new_cut_dim); + + CGAL_assertion(new_rd >= copy_rd); + dists[new_cut_dim] = diff1; + Node_with_distance *Upper_Child = + new Node_with_distance(node->upper(), new_rd, dists); + PriorityQueue.push(Upper_Child); + dists[new_cut_dim] = dst; + N=node->lower(); + + } + else { // compute new distance + new_rd=Orthogonal_distance_instance.new_distance(copy_rd,dst,diff2,new_cut_dim); + CGAL_assertion(new_rd >= copy_rd); + dists[new_cut_dim] = diff2; + Node_with_distance *Lower_Child = + new Node_with_distance(node->lower(), new_rd, dists); + PriorityQueue.push(Lower_Child); + dists[new_cut_dim] = dst; + N=node->upper(); + } + } + // n is a leaf + typename Tree::Leaf_node_const_handle node = + static_cast(N); + number_of_leaf_nodes_visited++; + if (node->size() > 0) { + for (typename Tree::iterator it=node->begin(); it != node->end(); it++) { + number_of_items_visited++; + FT distance_to_query_point= + Orthogonal_distance_instance.transformed_distance(query_point,*it); + Point_with_transformed_distance *NN_Candidate= + new Point_with_transformed_distance(*it,distance_to_query_point); + Item_PriorityQueue.push(NN_Candidate); + } + // old top of PriorityQueue has been processed, + // hence update rd + + if (!(PriorityQueue.empty())) { + rd = CGAL::cpp11::get<1>(*PriorityQueue.top()); + next_neighbour_found = + (multiplication_factor*rd > + Item_PriorityQueue.top()->second); + } + else // priority queue empty => last neighbour found + { + next_neighbour_found=true; + } + + number_of_neighbours_computed++; + } + } // next_neighbour_found or priority queue is empty + // in the latter case also the item priority quee is empty + } + + + void + Compute_the_next_furthest_neighbour() + { + // compute the next item + bool next_neighbour_found=false; + if (!(Item_PriorityQueue.empty())) { + next_neighbour_found= + (rd < multiplication_factor*Item_PriorityQueue.top()->second); + } + typename SearchTraits::Construct_cartesian_const_iterator_d construct_it=traits.construct_cartesian_const_iterator_d_object(); + typename SearchTraits::Cartesian_const_iterator_d query_point_it = construct_it(query_point); + // otherwise browse the tree further + while ((!next_neighbour_found) && (!PriorityQueue.empty())) { + Node_with_distance* The_node_top=PriorityQueue.top(); + Node_const_handle N= CGAL::cpp11::get<0>(*The_node_top); + dists = CGAL::cpp11::get<2>(*The_node_top); + PriorityQueue.pop(); + delete The_node_top; + FT copy_rd=rd; + while (!(N->is_leaf())) { // compute new distance + typename Tree::Internal_node_const_handle node = + static_cast(N); + number_of_internal_nodes_visited++; + int new_cut_dim=node->cutting_dimension(); + FT new_rd,dst = dists[new_cut_dim]; + FT val = *(query_point_it + new_cut_dim); + FT diff1 = val - node->upper_low_value(); + FT diff2 = val - node->lower_high_value(); + if (diff1 + diff2 < FT(0.0)) { + diff1 = val - node->upper_high_value(); + new_rd= + Orthogonal_distance_instance.new_distance(copy_rd,dst,diff1,new_cut_dim); + Node_with_distance *Lower_Child = + new Node_with_distance(node->lower(), copy_rd, dists); + PriorityQueue.push(Lower_Child); + N=node->upper(); + dists[new_cut_dim] = diff1; + copy_rd=new_rd; + + } + else { // compute new distance + diff2 = val - node->lower_low_value(); + new_rd=Orthogonal_distance_instance.new_distance(copy_rd,dst,diff2,new_cut_dim); + Node_with_distance *Upper_Child = + new Node_with_distance(node->upper(), copy_rd, dists); + PriorityQueue.push(Upper_Child); + N=node->lower(); + dists[new_cut_dim] = diff2; + copy_rd=new_rd; + } + } + // n is a leaf + typename Tree::Leaf_node_const_handle node = + static_cast(N); + number_of_leaf_nodes_visited++; + if (node->size() > 0) { + for (typename Tree::iterator it=node->begin(); it != node->end(); it++) { + number_of_items_visited++; + FT distance_to_query_point= + Orthogonal_distance_instance.transformed_distance(query_point,*it); + Point_with_transformed_distance *NN_Candidate= + new Point_with_transformed_distance(*it,distance_to_query_point); + Item_PriorityQueue.push(NN_Candidate); + } + // old top of PriorityQueue has been processed, + // hence update rd + + if (!(PriorityQueue.empty())) { + rd = CGAL::cpp11::get<1>(*PriorityQueue.top()); + next_neighbour_found = + (multiplication_factor*rd < + Item_PriorityQueue.top()->second); + } + else // priority queue empty => last neighbour found + { + next_neighbour_found=true; + } + + number_of_neighbours_computed++; + } + } // next_neighbour_found or priority queue is empty + // in the latter case also the item priority quee is empty + } + }; // class Iterator_implementaion + + + + + + + + + + public: + class iterator; + typedef iterator const_iterator; + + // constructor + Orthogonal_incremental_neighbor_search(const Tree& tree, + const Query_item& q, FT Eps = FT(0.0), + bool search_nearest=true, const Distance& tr=Distance()) + : m_tree(tree),m_query(q),m_dist(tr),m_Eps(Eps),m_search_nearest(search_nearest) + {} + + iterator + begin() const + { + return iterator(m_tree,m_query,m_dist,m_Eps,m_search_nearest); + } + + iterator + end() const + { + return iterator(); + } + + std::ostream& + statistics(std::ostream& s) + { + begin()->statistics(s); + return s; + } + + + + + class iterator { + + public: + + typedef std::input_iterator_tag iterator_category; + typedef Point_with_transformed_distance value_type; + typedef Point_with_transformed_distance* pointer; + typedef const Point_with_transformed_distance& reference; + typedef std::size_t size_type; + typedef std::ptrdiff_t difference_type; + typedef int distance_type; + + //class Iterator_implementation; + Iterator_implementation *Ptr_implementation; + + + public: + + // default constructor + iterator() + : Ptr_implementation(0) + {} + + int + the_number_of_items_visited() + { + return Ptr_implementation->number_of_items_visited; + } + + // constructor + iterator(const Tree& tree,const Query_item& q, const Distance& tr=Distance(), FT eps=FT(0.0), + bool search_nearest=true) + : Ptr_implementation(new Iterator_implementation(tree, q, tr, eps, search_nearest)) + {} + + // copy constructor + iterator(const iterator& Iter) + { + Ptr_implementation = Iter.Ptr_implementation; + if (Ptr_implementation != 0) Ptr_implementation->reference_count++; + } + + iterator& operator=(const iterator& Iter) + { + if (Ptr_implementation != Iter.Ptr_implementation){ + if (Ptr_implementation != 0 && --(Ptr_implementation->reference_count)==0) { + delete Ptr_implementation; + } + Ptr_implementation = Iter.Ptr_implementation; + if (Ptr_implementation != 0) Ptr_implementation->reference_count++; + } + return *this; + } + + + const Point_with_transformed_distance& + operator* () const + { + return *(*Ptr_implementation); + } + + // -> operator + const Point_with_transformed_distance* + operator-> () const + { + return &*(*Ptr_implementation); + } + + // prefix operator + iterator& + operator++() + { + ++(*Ptr_implementation); + return *this; + } + + // postfix operator + Object_wrapper + operator++(int) + { + return (*Ptr_implementation)++; + } + + + bool + operator==(const iterator& It) const + { + if ( + ((Ptr_implementation == 0) || + Ptr_implementation->Item_PriorityQueue.empty()) && + ((It.Ptr_implementation == 0) || + It.Ptr_implementation->Item_PriorityQueue.empty()) + ) + return true; + // else + return (Ptr_implementation == It.Ptr_implementation); + } + + bool + operator!=(const iterator& It) const + { + return !(*this == It); + } + + std::ostream& + statistics (std::ostream& s) + { + Ptr_implementation->statistics(s); + return s; + } + + ~iterator() + { + if (Ptr_implementation != 0) { + Ptr_implementation->reference_count--; + if (Ptr_implementation->reference_count==0) { + delete Ptr_implementation; + Ptr_implementation = 0; + } + } + } + + + }; // class iterator + + //data members + const Tree& m_tree; + Query_item m_query; + Distance m_dist; + FT m_Eps; + bool m_search_nearest; + }; // class + + template + void swap (typename Orthogonal_incremental_neighbor_search::iterator& x, + typename Orthogonal_incremental_neighbor_search::iterator& y) + { + typename Orthogonal_incremental_neighbor_search::iterator::Iterator_implementation + *tmp = x.Ptr_implementation; + x.Ptr_implementation = y.Ptr_implementation; + y.Ptr_implementation = tmp; + } + +} // namespace CGAL + +#endif // CGAL_ORTHOGONAL_INCREMENTAL_NEIGHBOR_SEARCH_H diff --git a/src/common/include/gudhi_patches/Tangential_complex_CGAL_patches.txt b/src/common/include/gudhi_patches/Tangential_complex_CGAL_patches.txt new file mode 100644 index 00000000..5b9581a0 --- /dev/null +++ b/src/common/include/gudhi_patches/Tangential_complex_CGAL_patches.txt @@ -0,0 +1,82 @@ +CGAL/Regular_triangulation_traits_adapter.h +CGAL/Triangulation_ds_vertex.h +CGAL/Triangulation_data_structure.h +CGAL/transforming_pair_iterator.h +CGAL/NewKernel_d/static_int.h +CGAL/NewKernel_d/Cartesian_LA_functors.h +CGAL/NewKernel_d/Cartesian_change_FT.h +CGAL/NewKernel_d/Wrapper/Vector_d.h +CGAL/NewKernel_d/Wrapper/Hyperplane_d.h +CGAL/NewKernel_d/Wrapper/Ref_count_obj.h +CGAL/NewKernel_d/Wrapper/Cartesian_wrap.h +CGAL/NewKernel_d/Wrapper/Point_d.h +CGAL/NewKernel_d/Wrapper/Segment_d.h +CGAL/NewKernel_d/Wrapper/Weighted_point_d.h +CGAL/NewKernel_d/Wrapper/Sphere_d.h +CGAL/NewKernel_d/Cartesian_per_dimension.h +CGAL/NewKernel_d/Kernel_object_converter.h +CGAL/NewKernel_d/KernelD_converter.h +CGAL/NewKernel_d/Vector/sse2.h +CGAL/NewKernel_d/Vector/avx4.h +CGAL/NewKernel_d/Vector/determinant_of_vectors_small_dim_internal.h +CGAL/NewKernel_d/Vector/determinant_of_iterator_to_points_from_points.h +CGAL/NewKernel_d/Vector/determinant_of_points_from_vectors.h +CGAL/NewKernel_d/Vector/array.h +CGAL/NewKernel_d/Vector/determinant_of_iterator_to_points_from_iterator_to_vectors.h +CGAL/NewKernel_d/Vector/determinant_of_iterator_to_vectors_from_vectors.h +CGAL/NewKernel_d/Vector/determinant_of_vectors_small_dim.h +CGAL/NewKernel_d/Vector/vector.h +CGAL/NewKernel_d/Vector/v2int.h +CGAL/NewKernel_d/Vector/mix.h +CGAL/NewKernel_d/Cartesian_static_filters.h +CGAL/NewKernel_d/Cartesian_LA_base.h +CGAL/NewKernel_d/Lazy_cartesian.h +CGAL/NewKernel_d/Coaffine.h +CGAL/NewKernel_d/store_kernel.h +CGAL/NewKernel_d/Dimension_base.h +CGAL/NewKernel_d/Kernel_3_interface.h +CGAL/NewKernel_d/Cartesian_complete.h +CGAL/NewKernel_d/Cartesian_base.h +CGAL/NewKernel_d/Cartesian_filter_K.h +CGAL/NewKernel_d/functor_tags.h +CGAL/NewKernel_d/Filtered_predicate2.h +CGAL/NewKernel_d/functor_properties.h +CGAL/NewKernel_d/Define_kernel_types.h +CGAL/NewKernel_d/LA_eigen/LA.h +CGAL/NewKernel_d/LA_eigen/constructors.h +CGAL/NewKernel_d/Types/Aff_transformation.h +CGAL/NewKernel_d/Types/Sphere.h +CGAL/NewKernel_d/Types/Hyperplane.h +CGAL/NewKernel_d/Types/Line.h +CGAL/NewKernel_d/Types/Ray.h +CGAL/NewKernel_d/Types/Iso_box.h +CGAL/NewKernel_d/Types/Weighted_point.h +CGAL/NewKernel_d/Types/Segment.h +CGAL/NewKernel_d/Kernel_d_interface.h +CGAL/NewKernel_d/utils.h +CGAL/NewKernel_d/Kernel_2_interface.h +CGAL/NewKernel_d/Cartesian_filter_NT.h +CGAL/NewKernel_d/function_objects_cartesian.h +CGAL/Convex_hull.h +CGAL/Triangulation_ds_full_cell.h +CGAL/Regular_triangulation.h +CGAL/Epick_d.h +CGAL/transforming_iterator.h +CGAL/iterator_from_indices.h +CGAL/Delaunay_triangulation.h +CGAL/IO/Triangulation_off_ostream.h +CGAL/typeset.h +CGAL/Triangulation_full_cell.h +CGAL/Triangulation.h +CGAL/internal/Static_or_dynamic_array.h +CGAL/internal/Combination_enumerator.h +CGAL/internal/Triangulation/utilities.h +CGAL/internal/Triangulation/Triangulation_ds_iterators.h +CGAL/internal/Triangulation/Dummy_TDS.h +CGAL/argument_swaps.h +CGAL/Epeck_d.h +CGAL/determinant_of_vectors.h +CGAL/TDS_full_cell_default_storage_policy.h +CGAL/TDS_full_cell_mirror_storage_policy.h +CGAL/Triangulation_face.h +CGAL/Triangulation_vertex.h -- cgit v1.2.3