diff options
author | vrouvrea <vrouvrea@636b058d-ea47-450e-bf9e-a15bfbe3eedb> | 2017-01-24 09:34:04 +0000 |
---|---|---|
committer | vrouvrea <vrouvrea@636b058d-ea47-450e-bf9e-a15bfbe3eedb> | 2017-01-24 09:34:04 +0000 |
commit | be37aaeadc31ed10ede53393237d939d4aa47c82 (patch) | |
tree | 6b990088ce68cc19b3db1b8789eab831e62dc747 /src/Persistent_cohomology | |
parent | 30a66f9431d059eed620e0535583e914fa5a3c74 (diff) | |
parent | cd9d7681e06f29a25cf5775a384ac2e07e34abe9 (diff) |
Merge last trunk modifications
git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/ST_cythonize@1989 636b058d-ea47-450e-bf9e-a15bfbe3eedb
Former-commit-id: 0743b9beb802b839357ecce17e11af5d4ef2a163
Diffstat (limited to 'src/Persistent_cohomology')
16 files changed, 546 insertions, 175 deletions
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/example/performance_rips_persistence.cpp b/src/Persistent_cohomology/benchmark/performance_rips_persistence.cpp index b4d282ac..ba752999 100644 --- a/src/Persistent_cohomology/example/performance_rips_persistence.cpp +++ b/src/Persistent_cohomology/benchmark/performance_rips_persistence.cpp @@ -20,20 +20,26 @@ * along with this program. If not, see <http://www.gnu.org/licenses/>. */ -#include <gudhi/reader_utils.h> -#include <gudhi/graph_simplicial_complex.h> +#include <gudhi/Rips_complex.h> #include <gudhi/distance_functions.h> #include <gudhi/Simplex_tree.h> #include <gudhi/Persistent_cohomology.h> #include <gudhi/Persistent_cohomology/Multi_field.h> #include <gudhi/Hasse_complex.h> +#include <gudhi/Points_off_io.h> #include <chrono> #include <string> #include <vector> -using namespace Gudhi; -using namespace Gudhi::persistent_cohomology; +// Types definition +using Simplex_tree = Gudhi::Simplex_tree<Gudhi::Simplex_tree_options_fast_persistence>; +using Filtration_value = Simplex_tree::Filtration_value; +using Rips_complex = Gudhi::rips_complex::Rips_complex<Filtration_value>; +using Field_Zp = Gudhi::persistent_cohomology::Field_Zp; +using Multi_field = Gudhi::persistent_cohomology::Multi_field; +using Point = std::vector<double>; +using Points_off_reader = Gudhi::Points_off_reader<Point>; /* Compute the persistent homology of the complex cpx with coefficients in Z/pZ. */ template< typename FilteredComplex> @@ -66,33 +72,29 @@ int main(int argc, char * argv[]) { int elapsed_sec; { - std::string filepoints = "../../../data/points/Kl.txt"; + 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 filepoints - typedef std::vector<double> Point_t; - std::vector< Point_t > points; - read_points(filepoints, points); + // 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(); - Graph_t prox_graph = compute_proximity_graph(points, threshold - , euclidean_distance<Point_t>); + 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<std::chrono::milliseconds>(end - start).count(); std::cout << "Compute Rips graph in " << elapsed_sec << " ms.\n"; // Construct the Rips complex in a Simplex Tree - Simplex_tree<Simplex_tree_options_fast_persistence> st; + 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); + rips_complex_from_file.create_complex(st, dim_max); end = std::chrono::system_clock::now(); elapsed_sec = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count(); @@ -120,7 +122,7 @@ int main(int argc, char * argv[]) { // Convert the simplex tree into a hasse diagram start = std::chrono::system_clock::now(); - Hasse_complex<> hcpx(st); + Gudhi::Hasse_complex<> hcpx(st); end = std::chrono::system_clock::now(); elapsed_sec = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count(); std::cout << "Convert the simplex tree into a Hasse diagram in " << elapsed_sec << " ms.\n"; @@ -152,7 +154,7 @@ timing_persistence(FilteredComplex & cpx int elapsed_sec; { start = std::chrono::system_clock::now(); - Persistent_cohomology< FilteredComplex, Field_Zp > pcoh(cpx); + Gudhi::persistent_cohomology::Persistent_cohomology< FilteredComplex, Field_Zp > pcoh(cpx); end = std::chrono::system_clock::now(); elapsed_sec = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count(); std::cout << " Initialize pcoh in " << elapsed_sec << " ms.\n"; @@ -186,7 +188,7 @@ timing_persistence(FilteredComplex & cpx int elapsed_sec; { start = std::chrono::system_clock::now(); - Persistent_cohomology< FilteredComplex, Multi_field > pcoh(cpx); + Gudhi::persistent_cohomology::Persistent_cohomology< FilteredComplex, Multi_field > pcoh(cpx); end = std::chrono::system_clock::now(); elapsed_sec = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count(); std::cout << " Initialize pcoh in " << elapsed_sec << " ms.\n"; diff --git a/src/Persistent_cohomology/doc/Intro_persistent_cohomology.h b/src/Persistent_cohomology/doc/Intro_persistent_cohomology.h index 433cfd3e..40dd3f93 100644 --- a/src/Persistent_cohomology/doc/Intro_persistent_cohomology.h +++ b/src/Persistent_cohomology/doc/Intro_persistent_cohomology.h @@ -144,17 +144,23 @@ namespace persistent_cohomology { We provide several example files: run these examples with -h for details on their use, and read the README file. \li <a href="_persistent_cohomology_2rips_persistence_8cpp-example.html"> -Persistent_cohomology/rips_persistence.cpp</a> computes the Rips complex of a point cloud and its persistence diagram. +Persistent_cohomology/rips_persistence.cpp</a> computes the Rips complex of a point cloud and outputs its persistence +diagram. +\code $> ./rips_persistence ../../data/points/tore3D_1307.off -r 0.25 -m 0.5 -d 3 -p 3 \endcode +\code The complex contains 177838 simplices + and has dimension 3 +3 0 0 inf +3 1 0.0983494 inf +3 1 0.104347 inf +3 2 0.138335 inf \endcode \li <a href="_persistent_cohomology_2rips_multifield_persistence_8cpp-example.html"> -Persistent_cohomology/rips_multifield_persistence.cpp</a> computes the Rips complex of a point cloud and its +Persistent_cohomology/rips_multifield_persistence.cpp</a> computes the Rips complex of a point cloud and outputs its persistence diagram with a family of field coefficients. -\li <a href="_persistent_cohomology_2performance_rips_persistence_8cpp-example.html"> -Persistent_cohomology/performance_rips_persistence.cpp</a> provides timings for the construction of the Rips complex -on a set of points sampling a Klein bottle in \f$\mathbb{R}^5\f$ with a simplex tree, its conversion to a -Hasse diagram and the computation of persistent homology and multi-field persistent homology for the -different representations. +\li <a href="_persistent_cohomology_2rips_distance_matrix_persistence_8cpp-example.html"> +Persistent_cohomology/rips_distance_matrix_persistence.cpp</a> computes the Rips complex of a distance matrix and +outputs its persistence diagram. \li <a href="_persistent_cohomology_2alpha_complex_3d_persistence_8cpp-example.html"> Persistent_cohomology/alpha_complex_3d_persistence.cpp</a> computes the persistent homology with diff --git a/src/Persistent_cohomology/example/CMakeLists.txt b/src/Persistent_cohomology/example/CMakeLists.txt index 758bd6b1..38d7e9a9 100644 --- a/src/Persistent_cohomology/example/CMakeLists.txt +++ b/src/Persistent_cohomology/example/CMakeLists.txt @@ -11,9 +11,15 @@ target_link_libraries(plain_homology ${Boost_SYSTEM_LIBRARY}) add_executable(persistence_from_simple_simplex_tree persistence_from_simple_simplex_tree.cpp) target_link_libraries(persistence_from_simple_simplex_tree ${Boost_SYSTEM_LIBRARY}) +add_executable(rips_distance_matrix_persistence rips_distance_matrix_persistence.cpp) +target_link_libraries(rips_distance_matrix_persistence ${Boost_SYSTEM_LIBRARY} ${Boost_PROGRAM_OPTIONS_LIBRARY}) + add_executable(rips_persistence rips_persistence.cpp) target_link_libraries(rips_persistence ${Boost_SYSTEM_LIBRARY} ${Boost_PROGRAM_OPTIONS_LIBRARY}) +add_executable(rips_persistence_step_by_step rips_persistence_step_by_step.cpp) +target_link_libraries(rips_persistence_step_by_step ${Boost_SYSTEM_LIBRARY} ${Boost_PROGRAM_OPTIONS_LIBRARY}) + add_executable(rips_persistence_via_boundary_matrix rips_persistence_via_boundary_matrix.cpp) target_link_libraries(rips_persistence_via_boundary_matrix ${Boost_SYSTEM_LIBRARY} ${Boost_PROGRAM_OPTIONS_LIBRARY}) @@ -23,15 +29,19 @@ target_link_libraries(persistence_from_file ${Boost_SYSTEM_LIBRARY} ${Boost_PROG if (TBB_FOUND) target_link_libraries(plain_homology ${TBB_LIBRARIES}) target_link_libraries(persistence_from_simple_simplex_tree ${TBB_LIBRARIES}) + target_link_libraries(rips_distance_matrix_persistence ${TBB_LIBRARIES}) target_link_libraries(rips_persistence ${TBB_LIBRARIES}) + target_link_libraries(rips_persistence_step_by_step ${TBB_LIBRARIES}) target_link_libraries(rips_persistence_via_boundary_matrix ${TBB_LIBRARIES}) target_link_libraries(persistence_from_file ${TBB_LIBRARIES}) endif() add_test(plain_homology ${CMAKE_CURRENT_BINARY_DIR}/plain_homology) add_test(persistence_from_simple_simplex_tree ${CMAKE_CURRENT_BINARY_DIR}/persistence_from_simple_simplex_tree 1 0) -add_test(rips_persistence_3 ${CMAKE_CURRENT_BINARY_DIR}/rips_persistence ${CMAKE_SOURCE_DIR}/data/points/Kl.txt -r 0.16 -d 3 -p 3 -m 100) -add_test(rips_persistence_via_boundary_matrix_3 ${CMAKE_CURRENT_BINARY_DIR}/rips_persistence_via_boundary_matrix ${CMAKE_SOURCE_DIR}/data/points/Kl.txt -r 0.16 -d 3 -p 3 -m 100) +add_test(rips_distance_matrix ${CMAKE_CURRENT_BINARY_DIR}/rips_distance_matrix_persistence ${CMAKE_SOURCE_DIR}/data/distance_matrix/full_square_distance_matrix.csv -r 1.0 -d 3 -p 3 -m 0) +add_test(rips_persistence_3 ${CMAKE_CURRENT_BINARY_DIR}/rips_persistence ${CMAKE_SOURCE_DIR}/data/points/tore3D_1307.off -r 0.25 -m 0.5 -d 3 -p 3) +add_test(rips_persistence_step_by_step_3 ${CMAKE_CURRENT_BINARY_DIR}/rips_persistence_step_by_step ${CMAKE_SOURCE_DIR}/data/points/tore3D_1307.off -r 0.25 -m 0.5 -d 3 -p 3) +add_test(rips_persistence_via_boundary_matrix_3 ${CMAKE_CURRENT_BINARY_DIR}/rips_persistence_via_boundary_matrix ${CMAKE_SOURCE_DIR}/data/points/Kl.off -r 0.16 -d 3 -p 3 -m 100) add_test(persistence_from_file_3_2_0 ${CMAKE_CURRENT_BINARY_DIR}/persistence_from_file ${CMAKE_SOURCE_DIR}/data/filtered_simplicial_complex/bunny_5000_complex.fsc -p 2 -m 0) add_test(persistence_from_file_3_3_100 ${CMAKE_CURRENT_BINARY_DIR}/persistence_from_file ${CMAKE_SOURCE_DIR}/data/filtered_simplicial_complex/bunny_5000_complex.fsc -p 3 -m 100) @@ -39,14 +49,10 @@ 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) + add_test(rips_multifield_persistence_2_71 ${CMAKE_CURRENT_BINARY_DIR}/rips_multifield_persistence ${CMAKE_SOURCE_DIR}/data/points/tore3D_1307.off -r 0.25 -m 0.5 -d 3 -p 2 -q 71) endif(GMPXX_FOUND) endif(GMP_FOUND) diff --git a/src/Persistent_cohomology/example/README b/src/Persistent_cohomology/example/README index 7803e5ab..2ac79398 100644 --- a/src/Persistent_cohomology/example/README +++ b/src/Persistent_cohomology/example/README @@ -10,13 +10,13 @@ Example of use of RIPS: Computation of the persistent homology with Z/2Z coefficients of the Rips complex on points sampling a Klein bottle: -./rips_persistence ../../data/points/Kl.txt -r 0.25 -d 3 -p 2 -m 100 +./rips_persistence ../../data/points/tore3D_1307.off -r 0.25 -m 0.5 -d 3 -p 2 output: -210 0 0 inf -210 1 0.0702103 inf -2 1 0.0702103 inf -2 2 0.159992 inf +2 0 0 inf +2 1 0.0983494 inf +2 1 0.104347 inf +2 2 0.138335 inf Every line is of this format: p1*...*pr dim b d @@ -29,31 +29,45 @@ where with Z/3Z coefficients: -./rips_persistence ../../data/points/Kl.txt -r 0.25 -d 3 -p 3 -m 100 +./rips_persistence ../../data/points/tore3D_1307.off -r 0.25 -m 0.5 -d 3 -p 3 output: -3 0 0 inf -3 1 0.0702103 inf +3 0 0 inf +3 1 0.0983494 inf +3 1 0.104347 inf +3 2 0.138335 inf and the computation with Z/2Z and Z/3Z coefficients simultaneously: -./rips_multifield_persistence ../../data/points/Kl.txt -r 0.25 -d 3 -p 2 -q 3 -m 100 +./rips_multifield_persistence ../../data/points/tore3D_1307.off -r 0.25 -m 0.12 -d 3 -p 2 -q 3 output: -6 0 0 inf -6 1 0.0702103 inf -2 1 0.0702103 inf -2 2 0.159992 inf +6 0 0 inf +6 1 0.0983494 inf +6 1 0.104347 inf +6 2 0.138335 inf +6 0 0 0.122545 +6 0 0 0.121171 +6 0 0 0.120964 +6 0 0 0.12057 +6 0 0 0.12047 +6 0 0 0.120414 and finally the computation with all Z/pZ for 2 <= p <= 71 (20 first prime numbers): - ./rips_multifield_persistence ../../data/points/Kl.txt -r 0.25 -d 3 -p 2 -q 71 -m 100 + ./rips_multifield_persistence ../../data/points/Kl.off -r 0.25 -m 0.5 -d 3 -p 2 -q 71 output: -557940830126698960967415390 0 0 inf -557940830126698960967415390 1 0.0702103 inf -2 1 0.0702103 inf -2 2 0.159992 inf +557940830126698960967415390 0 0 inf +557940830126698960967415390 1 0.0983494 inf +557940830126698960967415390 1 0.104347 inf +557940830126698960967415390 2 0.138335 inf +557940830126698960967415390 0 0 0.122545 +557940830126698960967415390 0 0 0.121171 +557940830126698960967415390 0 0 0.120964 +557940830126698960967415390 0 0 0.12057 +557940830126698960967415390 0 0 0.12047 +557940830126698960967415390 0 0 0.120414 *********************************************************************************************************************** Example of use of ALPHA: diff --git a/src/Persistent_cohomology/example/alpha_complex_3d_persistence.cpp b/src/Persistent_cohomology/example/alpha_complex_3d_persistence.cpp index 48fbb91a..978dc942 100644 --- a/src/Persistent_cohomology/example/alpha_complex_3d_persistence.cpp +++ b/src/Persistent_cohomology/example/alpha_complex_3d_persistence.cpp @@ -4,7 +4,7 @@ * * Author(s): Vincent Rouvreau * - * Copyright (C) 2014 INRIA Saclay (France) + * Copyright (C) 2014 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 @@ -64,6 +64,7 @@ typedef std::list<Alpha_shape_3::Vertex_handle> Vertex_list; // gudhi type definition typedef Gudhi::Simplex_tree<Gudhi::Simplex_tree_options_fast_persistence> ST; +typedef ST::Filtration_value Filtration_value; typedef ST::Vertex_handle Simplex_tree_vertex; typedef std::map<Alpha_shape_3::Vertex_handle, Simplex_tree_vertex > Alpha_shape_simplex_tree_map; typedef std::pair<Alpha_shape_3::Vertex_handle, Simplex_tree_vertex> 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 2412569a..9e84e91f 100644 --- a/src/Persistent_cohomology/example/alpha_complex_persistence.cpp +++ b/src/Persistent_cohomology/example/alpha_complex_persistence.cpp @@ -11,6 +11,9 @@ #include <string> #include <limits> // for numeric_limits +using Simplex_tree = Gudhi::Simplex_tree<>; +using Filtration_value = Simplex_tree::Filtration_value; + void program_options(int argc, char * argv[] , std::string & off_file_points , std::string & output_file_diag @@ -34,7 +37,7 @@ int main(int argc, char **argv) { using Kernel = CGAL::Epick_d< CGAL::Dynamic_dimension_tag >; Gudhi::alpha_complex::Alpha_complex<Kernel> alpha_complex_from_file(off_file_points); - Gudhi::Simplex_tree<> simplex; + Simplex_tree simplex; if (alpha_complex_from_file.create_complex(simplex, alpha_square_max_value)) { // ---------------------------------------------------------------------------- // Display information about the alpha complex @@ -48,7 +51,7 @@ int main(int argc, char **argv) { std::cout << "Simplex_tree dim: " << simplex.dimension() << std::endl; // Compute the persistence diagram of the complex - Gudhi::persistent_cohomology::Persistent_cohomology< Gudhi::Simplex_tree<>, + Gudhi::persistent_cohomology::Persistent_cohomology< Simplex_tree, Gudhi::persistent_cohomology::Field_Zp > pcoh(simplex); // initializes the coefficient field for homology pcoh.init_coefficients(coeff_field_characteristic); 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..dbc42706 100644 --- a/src/Persistent_cohomology/example/periodic_alpha_complex_3d_persistence.cpp +++ b/src/Persistent_cohomology/example/periodic_alpha_complex_3d_persistence.cpp @@ -4,7 +4,7 @@ * * Author(s): Vincent Rouvreau * - * Copyright (C) 2014 INRIA Saclay (France) + * Copyright (C) 2014 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 @@ -39,6 +39,7 @@ #include <utility> #include <list> #include <vector> +#include <cstdlib> // Traits using K = CGAL::Exact_predicates_inexact_constructions_kernel; @@ -70,6 +71,7 @@ using Vertex_list = std::list<Alpha_shape_3::Vertex_handle>; // gudhi type definition using ST = Gudhi::Simplex_tree<Gudhi::Simplex_tree_options_fast_persistence>; +using Filtration_value = ST::Filtration_value; using Simplex_tree_vertex = ST::Vertex_handle; using Alpha_shape_simplex_tree_map = std::map<Alpha_shape_3::Vertex_handle, Simplex_tree_vertex >; using Alpha_shape_simplex_tree_pair = std::pair<Alpha_shape_3::Vertex_handle, Simplex_tree_vertex>; @@ -136,19 +138,8 @@ int main(int argc, char * const argv[]) { usage(argv[0]); } - int coeff_field_characteristic = 0; - int returnedScanValue = sscanf(argv[3], "%d", &coeff_field_characteristic); - if ((returnedScanValue == EOF) || (coeff_field_characteristic <= 0)) { - std::cerr << "Error: " << argv[3] << " is not correct\n"; - usage(argv[0]); - } - - Filtration_value min_persistence = 0.0; - returnedScanValue = sscanf(argv[4], "%lf", &min_persistence); - if ((returnedScanValue == EOF) || (min_persistence < -1.0)) { - std::cerr << "Error: " << argv[4] << " is not correct\n"; - usage(argv[0]); - } + int coeff_field_characteristic = atoi(argv[3]); + Filtration_value min_persistence = strtof(argv[4], nullptr); // Read points from file std::string offInputFile(argv[1]); 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..7ca9410a 100644 --- a/src/Persistent_cohomology/example/persistence_from_simple_simplex_tree.cpp +++ b/src/Persistent_cohomology/example/persistence_from_simple_simplex_tree.cpp @@ -4,7 +4,7 @@ * * Author(s): Vincent Rouvreau * - * Copyright (C) 2014 INRIA Sophia Antipolis-Méditerranée (France) + * Copyright (C) 2014 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 @@ -29,13 +29,12 @@ #include <utility> #include <vector> -using namespace Gudhi; -using namespace Gudhi::persistent_cohomology; - -typedef std::vector< Vertex_handle > typeVectorVertex; -typedef std::pair<typeVectorVertex, Filtration_value> 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<Simplex_tree, Field_Zp >; +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_distance_matrix_persistence.cpp b/src/Persistent_cohomology/example/rips_distance_matrix_persistence.cpp new file mode 100644 index 00000000..8517e7f6 --- /dev/null +++ b/src/Persistent_cohomology/example/rips_distance_matrix_persistence.cpp @@ -0,0 +1,144 @@ +/* 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): Pawel Dlotko, 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 <http://www.gnu.org/licenses/>. + */ + +#include <gudhi/Rips_complex.h> +#include <gudhi/Simplex_tree.h> +#include <gudhi/Persistent_cohomology.h> +#include <gudhi/reader_utils.h> + +#include <boost/program_options.hpp> + +#include <string> +#include <vector> +#include <limits> // infinity + +// Types definition +using Simplex_tree = Gudhi::Simplex_tree<Gudhi::Simplex_tree_options_fast_persistence>; +using Filtration_value = Simplex_tree::Filtration_value; +using Rips_complex = Gudhi::rips_complex::Rips_complex<Filtration_value>; +using Field_Zp = Gudhi::persistent_cohomology::Field_Zp; +using Persistent_cohomology = Gudhi::persistent_cohomology::Persistent_cohomology<Simplex_tree, Field_Zp >; +using Distance_matrix = std::vector<std::vector<Filtration_value>>; + +void program_options(int argc, char * argv[] + , std::string & csv_matrix_file + , std::string & filediag + , Filtration_value & threshold + , int & dim_max + , int & p + , Filtration_value & min_persistence); + +int main(int argc, char * argv[]) { + std::string csv_matrix_file; + std::string filediag; + Filtration_value threshold; + int dim_max; + int p; + Filtration_value min_persistence; + + program_options(argc, argv, csv_matrix_file, filediag, threshold, dim_max, p, min_persistence); + + Distance_matrix distances = read_lower_triangular_matrix_from_csv_file<Filtration_value>(csv_matrix_file); + Rips_complex rips_complex_from_file(distances, threshold); + + // Construct the Rips complex in a Simplex Tree + Simplex_tree simplex_tree; + + 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(); + + // Compute the persistence diagram of the complex + 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 & csv_matrix_file + , std::string & filediag + , Filtration_value & threshold + , int & dim_max + , int & p + , Filtration_value & min_persistence) { + namespace po = boost::program_options; + po::options_description hidden("Hidden options"); + hidden.add_options() + ("input-file", po::value<std::string>(&csv_matrix_file), + "Name of file containing a distance matrix. Can be square or lower triangular matrix. Separator is ';'."); + + po::options_description visible("Allowed options", 100); + visible.add_options() + ("help,h", "produce help message") + ("output-file,o", po::value<std::string>(&filediag)->default_value(std::string()), + "Name of file in which the persistence diagram is written. Default print in std::cout") + ("max-edge-length,r", + po::value<Filtration_value>(&threshold)->default_value(std::numeric_limits<Filtration_value>::infinity()), + "Maximal length of an edge for the Rips complex construction.") + ("cpx-dimension,d", po::value<int>(&dim_max)->default_value(1), + "Maximal dimension of the Rips complex we want to compute.") + ("field-charac,p", po::value<int>(&p)->default_value(11), + "Characteristic p of the coefficient field Z/pZ for computing homology.") + ("min-persistence,m", po::value<Filtration_value>(&min_persistence), + "Minimal lifetime of homology feature to be recorded. Default is 0. Enter a negative value to see zero length intervals"); + + po::positional_options_description pos; + pos.add("input-file", 1); + + po::options_description all; + all.add(visible).add(hidden); + + po::variables_map vm; + po::store(po::command_line_parser(argc, argv). + options(all).positional(pos).run(), vm); + po::notify(vm); + + if (vm.count("help") || !vm.count("input-file")) { + std::cout << std::endl; + std::cout << "Compute the persistent homology with coefficient field Z/pZ \n"; + std::cout << "of a Rips complex defined on a set of distance matrix.\n \n"; + std::cout << "The output diagram contains one bar per line, written with the convention: \n"; + std::cout << " p dim b d \n"; + std::cout << "where dim is the dimension of the homological feature,\n"; + std::cout << "b and d are respectively the birth and death of the feature and \n"; + std::cout << "p is the characteristic of the field Z/pZ used for homology coefficients." << std::endl << std::endl; + + std::cout << "Usage: " << argv[0] << " [options] input-file" << std::endl << std::endl; + std::cout << visible << std::endl; + std::abort(); + } +} diff --git a/src/Persistent_cohomology/example/rips_multifield_persistence.cpp b/src/Persistent_cohomology/example/rips_multifield_persistence.cpp index c5cd775d..7674b5a5 100644 --- a/src/Persistent_cohomology/example/rips_multifield_persistence.cpp +++ b/src/Persistent_cohomology/example/rips_multifield_persistence.cpp @@ -4,7 +4,7 @@ * * Author(s): Clément Maria * - * Copyright (C) 2014 INRIA Sophia Antipolis-Méditerranée (France) + * Copyright (C) 2014 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,26 +20,29 @@ * along with this program. If not, see <http://www.gnu.org/licenses/>. */ -#include <gudhi/reader_utils.h> -#include <gudhi/graph_simplicial_complex.h> +#include <gudhi/Rips_complex.h> #include <gudhi/distance_functions.h> #include <gudhi/Simplex_tree.h> #include <gudhi/Persistent_cohomology.h> #include <gudhi/Persistent_cohomology/Multi_field.h> +#include <gudhi/Points_off_io.h> #include <boost/program_options.hpp> #include <string> #include <vector> -using namespace Gudhi; -using namespace Gudhi::persistent_cohomology; - -typedef int Vertex_handle; -typedef double Filtration_value; +// Types definition +using Simplex_tree = Gudhi::Simplex_tree<Gudhi::Simplex_tree_options_fast_persistence>; +using Filtration_value = Simplex_tree::Filtration_value; +using Rips_complex = Gudhi::rips_complex::Rips_complex<Filtration_value>; +using Multi_field = Gudhi::persistent_cohomology::Multi_field; +using Persistent_cohomology = Gudhi::persistent_cohomology::Persistent_cohomology<Simplex_tree, Multi_field >; +using Point = std::vector<double>; +using Points_off_reader = Gudhi::Points_off_reader<Point>; 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 +51,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,33 +59,26 @@ 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); - - // Extract the points from the file filepoints - typedef std::vector<double> Point_t; - std::vector< Point_t > points; - read_points(filepoints, points); + program_options(argc, argv, off_file_points, filediag, threshold, dim_max, min_p, max_p, min_persistence); - // Compute the proximity graph of the points - Graph_t prox_graph = compute_proximity_graph(points, threshold - , euclidean_distance<Point_t>); + 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 - typedef Simplex_tree<Simplex_tree_options_fast_persistence> 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); + Simplex_tree simplex_tree; + + 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 - st.initialize_filtration(); + simplex_tree.initialize_filtration(); // Compute the persistence diagram of the complex - Persistent_cohomology<ST, Multi_field > pcoh(st); + Persistent_cohomology pcoh(simplex_tree); // 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 @@ -98,7 +94,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 @@ -108,8 +104,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<std::string>(&filepoints), - "Name of file containing a point set. Format is one point per line: X1 ... Xd \n"); + ("input-file", po::value<std::string>(&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..c6378de7 100644 --- a/src/Persistent_cohomology/example/rips_persistence.cpp +++ b/src/Persistent_cohomology/example/rips_persistence.cpp @@ -4,7 +4,7 @@ * * Author(s): Clément Maria * - * Copyright (C) 2014 INRIA Sophia Antipolis-Méditerranée (France) + * Copyright (C) 2014 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,11 +20,11 @@ * along with this program. If not, see <http://www.gnu.org/licenses/>. */ -#include <gudhi/reader_utils.h> -#include <gudhi/graph_simplicial_complex.h> +#include <gudhi/Rips_complex.h> #include <gudhi/distance_functions.h> #include <gudhi/Simplex_tree.h> #include <gudhi/Persistent_cohomology.h> +#include <gudhi/Points_off_io.h> #include <boost/program_options.hpp> @@ -32,14 +32,17 @@ #include <vector> #include <limits> // infinity -using namespace Gudhi; -using namespace Gudhi::persistent_cohomology; - -typedef int Vertex_handle; -typedef double Filtration_value; +// Types definition +using Simplex_tree = Gudhi::Simplex_tree<Gudhi::Simplex_tree_options_fast_persistence>; +using Filtration_value = Simplex_tree::Filtration_value; +using Rips_complex = Gudhi::rips_complex::Rips_complex<Filtration_value>; +using Field_Zp = Gudhi::persistent_cohomology::Field_Zp; +using Persistent_cohomology = Gudhi::persistent_cohomology::Persistent_cohomology<Simplex_tree, Field_Zp >; +using Point = std::vector<double>; +using Points_off_reader = Gudhi::Points_off_reader<Point>; 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,40 +50,30 @@ 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<double> 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<Point_t>); + 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 - typedef Simplex_tree<Simplex_tree_options_fast_persistence> 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); + Simplex_tree simplex_tree; - std::cout << "The complex contains " << st.num_simplices() << " simplices \n"; - std::cout << " and has dimension " << st.dimension() << " \n"; + 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 - st.initialize_filtration(); + simplex_tree.initialize_filtration(); // Compute the persistence diagram of the complex - persistent_cohomology::Persistent_cohomology<ST, Field_Zp > pcoh(st); + Persistent_cohomology pcoh(simplex_tree); // initializes the coefficient field for homology pcoh.init_coefficients(p); @@ -99,7 +92,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 @@ -108,15 +101,16 @@ 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<std::string>(&filepoints), - "Name of file containing a point set. Format is one point per line: X1 ... Xd "); + ("input-file", po::value<std::string>(&off_file_points), + "Name of an OFF file containing a point set.\n"); po::options_description visible("Allowed options", 100); visible.add_options() ("help,h", "produce help message") ("output-file,o", po::value<std::string>(&filediag)->default_value(std::string()), "Name of file in which the persistence diagram is written. Default print in std::cout") - ("max-edge-length,r", po::value<Filtration_value>(&threshold)->default_value(std::numeric_limits<Filtration_value>::infinity()), + ("max-edge-length,r", + po::value<Filtration_value>(&threshold)->default_value(std::numeric_limits<Filtration_value>::infinity()), "Maximal length of an edge for the Rips complex construction.") ("cpx-dimension,d", po::value<int>(&dim_max)->default_value(1), "Maximal dimension of the Rips complex we want to compute.") diff --git a/src/Persistent_cohomology/example/rips_persistence_step_by_step.cpp b/src/Persistent_cohomology/example/rips_persistence_step_by_step.cpp new file mode 100644 index 00000000..c8f0921a --- /dev/null +++ b/src/Persistent_cohomology/example/rips_persistence_step_by_step.cpp @@ -0,0 +1,210 @@ +/* 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 <http://www.gnu.org/licenses/>. + */ + +#include <gudhi/graph_simplicial_complex.h> +#include <gudhi/distance_functions.h> +#include <gudhi/Simplex_tree.h> +#include <gudhi/Persistent_cohomology.h> +#include <gudhi/Points_off_io.h> + +#include <boost/program_options.hpp> + +#include <string> +#include <vector> +#include <limits> // infinity +#include <utility> // for pair +#include <map> + +// Types definition +using Simplex_tree = Gudhi::Simplex_tree<Gudhi::Simplex_tree_options_fast_persistence>; +using Vertex_handle = Simplex_tree::Vertex_handle; +using Filtration_value = Simplex_tree::Filtration_value; +using Graph_t = boost::adjacency_list < boost::vecS, boost::vecS, boost::undirectedS +, boost::property < vertex_filtration_t, Filtration_value > +, boost::property < edge_filtration_t, Filtration_value > +>; +using Edge_t = std::pair< Vertex_handle, Vertex_handle >; + +template< typename InputPointRange, typename Distance > +Graph_t compute_proximity_graph(InputPointRange &points, Filtration_value threshold, Distance distance); + +using Field_Zp = Gudhi::persistent_cohomology::Field_Zp; +using Persistent_cohomology = Gudhi::persistent_cohomology::Persistent_cohomology<Simplex_tree, Field_Zp >; +using Point = std::vector<double>; +using Points_off_reader = Gudhi::Points_off_reader<Point>; + +void program_options(int argc, char * argv[] + , std::string & off_file_points + , std::string & filediag + , Filtration_value & threshold + , int & dim_max + , int & p + , Filtration_value & min_persistence); + +int main(int argc, char * argv[]) { + std::string off_file_points; + std::string filediag; + Filtration_value threshold; + int dim_max; + int p; + Filtration_value min_persistence; + + program_options(argc, argv, off_file_points, filediag, threshold, dim_max, p, min_persistence); + + // Extract the points from the file filepoints + Points_off_reader off_reader(off_file_points); + + // Compute the proximity graph of the points + Graph_t prox_graph = compute_proximity_graph(off_reader.get_point_cloud(), threshold + , Euclidean_distance()); + + // Construct the Rips complex in a Simplex Tree + Simplex_tree 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 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(); + } + + return 0; +} + +void program_options(int argc, char * argv[] + , std::string & off_file_points + , std::string & filediag + , Filtration_value & threshold + , int & dim_max + , int & p + , Filtration_value & min_persistence) { + namespace po = boost::program_options; + po::options_description hidden("Hidden options"); + hidden.add_options() + ("input-file", po::value<std::string>(&off_file_points), + "Name of an OFF file containing a point set.\n"); + + po::options_description visible("Allowed options", 100); + visible.add_options() + ("help,h", "produce help message") + ("output-file,o", po::value<std::string>(&filediag)->default_value(std::string()), + "Name of file in which the persistence diagram is written. Default print in std::cout") + ("max-edge-length,r", + po::value<Filtration_value>(&threshold)->default_value(std::numeric_limits<Filtration_value>::infinity()), + "Maximal length of an edge for the Rips complex construction.") + ("cpx-dimension,d", po::value<int>(&dim_max)->default_value(1), + "Maximal dimension of the Rips complex we want to compute.") + ("field-charac,p", po::value<int>(&p)->default_value(11), + "Characteristic p of the coefficient field Z/pZ for computing homology.") + ("min-persistence,m", po::value<Filtration_value>(&min_persistence), + "Minimal lifetime of homology feature to be recorded. Default is 0. Enter a negative value to see zero length intervals"); + + po::positional_options_description pos; + pos.add("input-file", 1); + + po::options_description all; + all.add(visible).add(hidden); + + po::variables_map vm; + po::store(po::command_line_parser(argc, argv). + options(all).positional(pos).run(), vm); + po::notify(vm); + + if (vm.count("help") || !vm.count("input-file")) { + std::cout << std::endl; + std::cout << "Compute the persistent homology with coefficient field Z/pZ \n"; + std::cout << "of a Rips complex defined on a set of input points.\n \n"; + std::cout << "The output diagram contains one bar per line, written with the convention: \n"; + std::cout << " p dim b d \n"; + std::cout << "where dim is the dimension of the homological feature,\n"; + std::cout << "b and d are respectively the birth and death of the feature and \n"; + std::cout << "p is the characteristic of the field Z/pZ used for homology coefficients." << std::endl << std::endl; + + std::cout << "Usage: " << argv[0] << " [options] input-file" << std::endl << std::endl; + std::cout << visible << std::endl; + std::abort(); + } +} + +/** Output 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. + * + * The type PointCloud furnishes .begin() and .end() methods, that return + * iterators with value_type Point. + */ +template< typename InputPointRange, typename Distance > +Graph_t compute_proximity_graph(InputPointRange &points, Filtration_value threshold, Distance distance) { + std::vector< Edge_t > edges; + std::vector< Filtration_value > edges_fil; + + Vertex_handle idx_u, idx_v; + Filtration_value fil; + idx_u = 0; + for (auto it_u = points.begin(); it_u != points.end(); ++it_u) { + idx_v = idx_u + 1; + for (auto it_v = it_u + 1; it_v != points.end(); ++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; + } + + Graph_t skel_graph(edges.begin() + , edges.end() + , edges_fil.begin() + , idx_u); // number of points labeled from 0 to idx_u-1 + + auto vertex_prop = boost::get(vertex_filtration_t(), skel_graph); + + boost::graph_traits<Graph_t>::vertex_iterator vi, vi_end; + for (std::tie(vi, vi_end) = boost::vertices(skel_graph); + vi != vi_end; ++vi) { + boost::put(vertex_prop, *vi, 0.); + } + + return skel_graph; +} 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..63da9847 100644 --- a/src/Persistent_cohomology/example/rips_persistence_via_boundary_matrix.cpp +++ b/src/Persistent_cohomology/example/rips_persistence_via_boundary_matrix.cpp @@ -4,8 +4,7 @@ * * Author(s): Clément Maria, Marc Glisse * - * Copyright (C) 2014 INRIA Sophia Antipolis-Méditerranée (France), - * 2015 INRIA Saclay Île de France) + * Copyright (C) 2014 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 @@ -21,12 +20,12 @@ * along with this program. If not, see <http://www.gnu.org/licenses/>. */ -#include <gudhi/reader_utils.h> -#include <gudhi/graph_simplicial_complex.h> -#include <gudhi/distance_functions.h> #include <gudhi/Simplex_tree.h> #include <gudhi/Persistent_cohomology.h> +#include <gudhi/Rips_complex.h> #include <gudhi/Hasse_complex.h> +#include <gudhi/Points_off_io.h> +#include <gudhi/distance_functions.h> #include <boost/program_options.hpp> @@ -44,14 +43,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<Filtration_value>; +using Field_Zp = Gudhi::persistent_cohomology::Field_Zp; +using Point = std::vector<double>; +using Points_off_reader = Gudhi::Points_off_reader<Point>; 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 +60,21 @@ 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<double> 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<Point_t>); + 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 +91,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 +101,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 +118,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 +127,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<std::string>(&filepoints), + ("input-file", po::value<std::string>(&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/include/gudhi/Persistent_cohomology.h b/src/Persistent_cohomology/include/gudhi/Persistent_cohomology.h index 8e7901cc..672fda48 100644 --- a/src/Persistent_cohomology/include/gudhi/Persistent_cohomology.h +++ b/src/Persistent_cohomology/include/gudhi/Persistent_cohomology.h @@ -300,8 +300,7 @@ class Persistent_cohomology { // with multiplicity. We used to sum the coefficients directly in // annotations_in_boundary by using a map, we now do it later. typedef std::pair<Column *, int> annotation_t; - // Danger: not thread-safe! - static std::vector<annotation_t> annotations_in_boundary; + thread_local std::vector<annotation_t> annotations_in_boundary; annotations_in_boundary.clear(); int sign = 1 - 2 * (dim_sigma % 2); // \in {-1,1} provides the sign in the // alternate sum in the boundary. 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<Filtration_value> (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 |