diff options
Diffstat (limited to 'src')
45 files changed, 1921 insertions, 235 deletions
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index e26b2d25..db9c40fc 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -112,6 +112,7 @@ else() add_subdirectory(example/Bitmap_cubical_complex) add_subdirectory(example/Witness_complex) add_subdirectory(example/Alpha_complex) + add_subdirectory(example/Rips_complex) add_subdirectory(example/Spatial_searching) add_subdirectory(example/Subsampling) add_subdirectory(example/Tangential_complex) diff --git a/src/Doxyfile b/src/Doxyfile index 943869ad..e7c1a0e1 100644 --- a/src/Doxyfile +++ b/src/Doxyfile @@ -846,6 +846,7 @@ IMAGE_PATH = doc/Skeleton_blocker/ \ doc/Persistent_cohomology/ \ doc/Witness_complex/ \ doc/Bitmap_cubical_complex/ \ + doc/Rips_complex/ \ doc/Subsampling/ \ doc/Spatial_searching/ \ doc/Tangential_complex/ diff --git a/src/GudhUI/utils/Persistence_compute.h b/src/GudhUI/utils/Persistence_compute.h index 97165490..f7048c28 100644 --- a/src/GudhUI/utils/Persistence_compute.h +++ b/src/GudhUI/utils/Persistence_compute.h @@ -29,6 +29,7 @@ #include <gudhi/Simplex_tree.h> #include <gudhi/distance_functions.h> #include <gudhi/Persistent_cohomology.h> +#include <gudhi/Rips_complex.h> #include <vector> @@ -69,22 +70,26 @@ template<typename SkBlComplex> 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<Filtration_value>; + using Field_Zp = Gudhi::persistent_cohomology::Field_Zp; + using Persistent_cohomology = Gudhi::persistent_cohomology::Persistent_cohomology<Simplex_tree, Field_Zp>; + + Rips_complex rips_complex(points, params.threshold, euclidean_distance<Filtration_value, Point_t>); - Graph_t prox_graph = compute_proximity_graph(points, params.threshold, euclidean_distance<Point_t>); - 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/example/performance_rips_persistence.cpp b/src/Persistent_cohomology/benchmark/performance_rips_persistence.cpp index b4d282ac..a9eab5dd 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,30 @@ 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<Filtration_value, Point>); 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 +123,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 +155,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 +189,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/example/CMakeLists.txt b/src/Persistent_cohomology/example/CMakeLists.txt index 758bd6b1..e50e8ca9 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<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..bca4b66c 100644 --- a/src/Persistent_cohomology/example/alpha_complex_persistence.cpp +++ b/src/Persistent_cohomology/example/alpha_complex_persistence.cpp @@ -11,6 +11,8 @@ #include <string> #include <limits> // 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/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<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>; 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 <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</*Gudhi::Simplex_tree_options_fast_persistence*/>; +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_multifield_persistence.cpp b/src/Persistent_cohomology/example/rips_multifield_persistence.cpp index c5cd775d..8e5c83a1 100644 --- a/src/Persistent_cohomology/example/rips_multifield_persistence.cpp +++ b/src/Persistent_cohomology/example/rips_multifield_persistence.cpp @@ -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,49 +59,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); - - // 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<Filtration_value, Point>); // 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); - - // Sort the simplices in the order of the filtration - st.initialize_filtration(); - - // Compute the persistence diagram of the complex - Persistent_cohomology<ST, Multi_field > 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(); + Simplex_tree simplex_tree; + + 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(); + + // Compute the persistence diagram of the complex + 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 +105,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..b74d0094 100644 --- a/src/Persistent_cohomology/example/rips_persistence.cpp +++ b/src/Persistent_cohomology/example/rips_persistence.cpp @@ -20,11 +20,12 @@ * along with this program. If not, see <http://www.gnu.org/licenses/>. */ -#include <gudhi/reader_utils.h> +#include <gudhi/Rips_complex.h> #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> @@ -32,14 +33,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,59 +51,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<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<Filtration_value, Point>); // 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); - - 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<ST, Field_Zp > 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(); + Simplex_tree simplex_tree; + + 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(); + + // 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 & filepoints + , std::string & off_file_points , std::string & filediag , Filtration_value & threshold , int & dim_max @@ -108,7 +103,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/example/rips_persistence_via_boundary_matrix.cpp b/src/Persistent_cohomology/example/rips_persistence_via_boundary_matrix.cpp index 4c6656f5..ee6577f4 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 <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 +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<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 +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<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<Filtration_value, Point>); // 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<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/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 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 <http://www.gnu.org/licenses/>. + */ + +#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<class OneSkeletonGraph> + 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..bd9f5fad --- /dev/null +++ b/src/Rips_complex/doc/Intro_rips_complex.h @@ -0,0 +1,105 @@ +/* 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 <http://www.gnu.org/licenses/>. + */ + +#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 and Vincent Rouvreau + * + * @{ + * + * \section ripsdefinition Definition + * + * Rips_complex + * <a target="_blank" href="https://en.wikipedia.org/wiki/Vietoris%E2%80%93Rips_complex">(Wikipedia)</a> is a + * <a target="_blank" href="https://en.wikipedia.org/wiki/Simplicial_complex">simplicial complex</a> + * 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_one_skeleton_rips_from_points.cpp + * + * When launching (rips maximal distance between 2 points is 12.0, is expanded until dimension 1 - one skeleton graph + * with other words): + * + * \code $> ./oneskeletonripspoints 12.0 + * \endcode + * + * the program output is: + * + * \include Rips_complex/one_skeleton_rips_points_for_doc.txt + * + * \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. + * 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/full_skeleton_rips_points_for_doc.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 @@ +<?xml version="1.0"?> +<!DOCTYPE ipe SYSTEM "ipe.dtd"> +<ipe version="70107" creator="Ipe 7.1.10"> +<info created="D:20150603143945" modified="D:20160928121844"/> +<ipestyle name="basic"> +<symbol name="arrow/arc(spx)"> +<path stroke="sym-stroke" fill="sym-stroke" pen="sym-pen"> +0 0 m +-1 0.333 l +-1 -0.333 l +h +</path> +</symbol> +<symbol name="arrow/farc(spx)"> +<path stroke="sym-stroke" fill="white" pen="sym-pen"> +0 0 m +-1 0.333 l +-1 -0.333 l +h +</path> +</symbol> +<symbol name="mark/circle(sx)" transformations="translations"> +<path fill="sym-stroke"> +0.6 0 0 0.6 0 0 e +0.4 0 0 0.4 0 0 e +</path> +</symbol> +<symbol name="mark/disk(sx)" transformations="translations"> +<path fill="sym-stroke"> +0.6 0 0 0.6 0 0 e +</path> +</symbol> +<symbol name="mark/fdisk(sfx)" transformations="translations"> +<group> +<path fill="sym-fill"> +0.5 0 0 0.5 0 0 e +</path> +<path fill="sym-stroke" fillrule="eofill"> +0.6 0 0 0.6 0 0 e +0.4 0 0 0.4 0 0 e +</path> +</group> +</symbol> +<symbol name="mark/box(sx)" transformations="translations"> +<path fill="sym-stroke" fillrule="eofill"> +-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 +</path> +</symbol> +<symbol name="mark/square(sx)" transformations="translations"> +<path fill="sym-stroke"> +-0.6 -0.6 m +0.6 -0.6 l +0.6 0.6 l +-0.6 0.6 l +h +</path> +</symbol> +<symbol name="mark/fsquare(sfx)" transformations="translations"> +<group> +<path fill="sym-fill"> +-0.5 -0.5 m +0.5 -0.5 l +0.5 0.5 l +-0.5 0.5 l +h +</path> +<path fill="sym-stroke" fillrule="eofill"> +-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 +</path> +</group> +</symbol> +<symbol name="mark/cross(sx)" transformations="translations"> +<group> +<path fill="sym-stroke"> +-0.43 -0.57 m +0.57 0.43 l +0.43 0.57 l +-0.57 -0.43 l +h +</path> +<path fill="sym-stroke"> +-0.43 0.57 m +0.57 -0.43 l +0.43 -0.57 l +-0.57 0.43 l +h +</path> +</group> +</symbol> +<symbol name="arrow/fnormal(spx)"> +<path stroke="sym-stroke" fill="white" pen="sym-pen"> +0 0 m +-1 0.333 l +-1 -0.333 l +h +</path> +</symbol> +<symbol name="arrow/pointed(spx)"> +<path stroke="sym-stroke" fill="sym-stroke" pen="sym-pen"> +0 0 m +-1 0.333 l +-0.8 0 l +-1 -0.333 l +h +</path> +</symbol> +<symbol name="arrow/fpointed(spx)"> +<path stroke="sym-stroke" fill="white" pen="sym-pen"> +0 0 m +-1 0.333 l +-0.8 0 l +-1 -0.333 l +h +</path> +</symbol> +<symbol name="arrow/linear(spx)"> +<path stroke="sym-stroke" pen="sym-pen"> +-1 0.333 m +0 0 l +-1 -0.333 l +</path> +</symbol> +<symbol name="arrow/fdouble(spx)"> +<path stroke="sym-stroke" fill="white" pen="sym-pen"> +0 0 m +-1 0.333 l +-1 -0.333 l +h +-1 0 m +-2 0.333 l +-2 -0.333 l +h +</path> +</symbol> +<symbol name="arrow/double(spx)"> +<path stroke="sym-stroke" fill="sym-stroke" pen="sym-pen"> +0 0 m +-1 0.333 l +-1 -0.333 l +h +-1 0 m +-2 0.333 l +-2 -0.333 l +h +</path> +</symbol> +<pen name="heavier" value="0.8"/> +<pen name="fat" value="1.2"/> +<pen name="ultrafat" value="2"/> +<symbolsize name="large" value="5"/> +<symbolsize name="small" value="2"/> +<symbolsize name="tiny" value="1.1"/> +<arrowsize name="large" value="10"/> +<arrowsize name="small" value="5"/> +<arrowsize name="tiny" value="3"/> +<color name="red" value="1 0 0"/> +<color name="green" value="0 1 0"/> +<color name="blue" value="0 0 1"/> +<color name="yellow" value="1 1 0"/> +<color name="orange" value="1 0.647 0"/> +<color name="gold" value="1 0.843 0"/> +<color name="purple" value="0.627 0.125 0.941"/> +<color name="gray" value="0.745"/> +<color name="brown" value="0.647 0.165 0.165"/> +<color name="navy" value="0 0 0.502"/> +<color name="pink" value="1 0.753 0.796"/> +<color name="seagreen" value="0.18 0.545 0.341"/> +<color name="turquoise" value="0.251 0.878 0.816"/> +<color name="violet" value="0.933 0.51 0.933"/> +<color name="darkblue" value="0 0 0.545"/> +<color name="darkcyan" value="0 0.545 0.545"/> +<color name="darkgray" value="0.663"/> +<color name="darkgreen" value="0 0.392 0"/> +<color name="darkmagenta" value="0.545 0 0.545"/> +<color name="darkorange" value="1 0.549 0"/> +<color name="darkred" value="0.545 0 0"/> +<color name="lightblue" value="0.678 0.847 0.902"/> +<color name="lightcyan" value="0.878 1 1"/> +<color name="lightgray" value="0.827"/> +<color name="lightgreen" value="0.565 0.933 0.565"/> +<color name="lightyellow" value="1 1 0.878"/> +<dashstyle name="dashed" value="[4] 0"/> +<dashstyle name="dotted" value="[1 3] 0"/> +<dashstyle name="dash dotted" value="[4 2 1 2] 0"/> +<dashstyle name="dash dot dotted" value="[4 2 1 2 1 2] 0"/> +<textsize name="large" value="\large"/> +<textsize name="small" value="\small"/> +<textsize name="tiny" value="\tiny"/> +<textsize name="Large" value="\Large"/> +<textsize name="LARGE" value="\LARGE"/> +<textsize name="huge" value="\huge"/> +<textsize name="Huge" value="\Huge"/> +<textsize name="footnote" value="\footnotesize"/> +<textstyle name="center" begin="\begin{center}" end="\end{center}"/> +<textstyle name="itemize" begin="\begin{itemize}" end="\end{itemize}"/> +<textstyle name="item" begin="\begin{itemize}\item{}" end="\end{itemize}"/> +<gridsize name="4 pts" value="4"/> +<gridsize name="8 pts (~3 mm)" value="8"/> +<gridsize name="16 pts (~6 mm)" value="16"/> +<gridsize name="32 pts (~12 mm)" value="32"/> +<gridsize name="10 pts (~3.5 mm)" value="10"/> +<gridsize name="20 pts (~7 mm)" value="20"/> +<gridsize name="14 pts (~5 mm)" value="14"/> +<gridsize name="28 pts (~10 mm)" value="28"/> +<gridsize name="56 pts (~20 mm)" value="56"/> +<anglesize name="90 deg" value="90"/> +<anglesize name="60 deg" value="60"/> +<anglesize name="45 deg" value="45"/> +<anglesize name="30 deg" value="30"/> +<anglesize name="22.5 deg" value="22.5"/> +<tiling name="falling" angle="-60" step="4" width="1"/> +<tiling name="rising" angle="30" step="4" width="1"/> +</ipestyle> +<page> +<layer name="alpha"/> +<view layers="alpha" active="alpha"/> +<path layer="alpha" matrix="1 0 0 1 0 -8" fill="darkblue"> +109.771 601.912 m +159.595 601.797 l +140.058 541.915 l +h +</path> +<path matrix="1 0 0 1 0 -8" fill="darkblue"> +79.8776 552.169 m +109.756 601.699 l +139.812 542.209 l +h +</path> +<path matrix="1 0 0 1 0 -8" fill="lightblue"> +69.8453 682.419 m +159.925 712.208 l +90.12 732.039 l +h +</path> +<text matrix="1 0 0 1 -230.178 14.1775" transformations="translations" pos="380 530" stroke="seagreen" type="label" width="68.836" height="8.307" depth="2.32" valign="baseline" size="large">Rips complex</text> +<text matrix="1 0 0 1 -212.333 10.6762" transformations="translations" pos="282.952 524.893" stroke="black" type="label" width="4.981" height="6.42" depth="0" valign="baseline">0</text> +<text matrix="1 0 0 1 -210.178 14.1775" transformations="translations" pos="352.708 510.349" stroke="black" type="label" width="4.981" height="6.42" depth="0" valign="baseline">1</text> +<text matrix="1 0 0 1 -210.178 14.1775" transformations="translations" pos="310.693 578.759" stroke="black" type="label" width="4.981" height="6.42" depth="0" valign="baseline">2</text> +<text matrix="1 0 0 1 -210.178 14.1775" transformations="translations" pos="375.332 578.49" stroke="black" type="label" width="4.981" height="6.42" depth="0" valign="baseline">3</text> +<text matrix="1 0 0 1 -210.178 14.1775" transformations="translations" pos="272.179 660.635" stroke="black" type="label" width="4.981" height="6.42" depth="0" valign="baseline">4</text> +<text matrix="1 0 0 1 -209.478 4.0238" transformations="translations" pos="296.419 724.197" stroke="black" type="label" width="4.981" height="6.42" depth="0" valign="baseline">5</text> +<text matrix="1 0 0 1 -210.178 14.1775" transformations="translations" pos="375.332 689.453" stroke="black" type="label" width="4.981" height="6.42" depth="0" valign="baseline">6</text> +<path matrix="1 0 0 1 29.8225 14.1775" stroke="black" pen="heavier"> +60 710 m +40 660 l +</path> +<path matrix="1 0 0 1 29.8225 14.1775" stroke="black" pen="heavier"> +40 660 m +130 690 l +</path> +<path matrix="1 0 0 1 29.8225 14.1775" stroke="black" pen="heavier"> +130 690 m +60 710 l +</path> +<path matrix="1 0 0 1 29.8225 14.1775" stroke="black" pen="heavier"> +40 660 m +80 580 l +</path> +<path matrix="1 0 0 1 29.8225 14.1775" stroke="black" pen="heavier"> +80 580 m +130 580 l +130 580 l +</path> +<path matrix="1 0 0 1 29.8225 14.1775" stroke="black" pen="heavier"> +130 580 m +110 520 l +</path> +<path matrix="1 0 0 1 29.8225 14.1775" stroke="black" pen="heavier"> +110 520 m +50 530 l +50 530 l +50 530 l +</path> +<path matrix="1 0 0 1 29.8225 14.1775" stroke="black" pen="heavier"> +50 530 m +80 580 l +</path> +<path matrix="1 0 0 1 29.8225 14.1775" stroke="black" pen="heavier"> +130 580 m +130 690 l +</path> +<use matrix="1 0 0 1 -209.478 4.0238" name="mark/fdisk(sfx)" pos="300 720" size="normal" stroke="black" fill="white"/> +<use matrix="1 0 0 1 -210.178 14.1775" name="mark/fdisk(sfx)" pos="280 660" size="normal" stroke="black" fill="white"/> +<use matrix="1 0 0 1 -210.178 14.1775" name="mark/fdisk(sfx)" pos="370 690" size="normal" stroke="black" fill="white"/> +<use matrix="1 0 0 1 -210.178 14.1775" name="mark/fdisk(sfx)" pos="370 580" size="normal" stroke="black" fill="white"/> +<use matrix="1 0 0 1 -210.178 14.1775" name="mark/fdisk(sfx)" pos="290 530" size="normal" stroke="black" fill="white"/> +<path matrix="1 0 0 1 -40 -16" stroke="black" pen="heavier"> +150.038 609.9 m +179.929 549.727 l +</path> +<use matrix="1 0 0 1 -210.178 14.1775" name="mark/fdisk(sfx)" pos="320 580" size="normal" stroke="black" fill="white"/> +<use matrix="1 0 0 1 -210.178 14.1775" name="mark/fdisk(sfx)" pos="350 520" size="normal" stroke="black" fill="white"/> +<path stroke="black" pen="heavier"> +158.7 593.269 m +81.4925 544.805 l +</path> +<path matrix="1 0 0 1 -17.9662 -17.9662" stroke="gray"> +256.324 639.958 m +370.055 639.958 l +</path> +<path matrix="1 0 0 1 -17.9662 -17.9662" stroke="gray"> +56.8567 0 0 56.8567 313.217 639.756 e +</path> +<use matrix="1 0 0 1 52.1387 -98.0941" name="mark/fdisk(sfx)" pos="300 720" size="normal" stroke="gray" fill="white"/> +<use matrix="1 0 0 1 -61.4926 -98.0942" name="mark/fdisk(sfx)" pos="300 720" size="normal" stroke="gray" fill="white"/> +<text matrix="1 0 0 1 -26.6167 -33.2708" transformations="translations" pos="295.735 657.944" stroke="gray" type="label" width="63.374" height="6.926" depth="1.93" valign="baseline">Rips threshold</text> +</page> +</ipe> diff --git a/src/Rips_complex/doc/rips_complex_representation.png b/src/Rips_complex/doc/rips_complex_representation.png Binary files differnew file mode 100644 index 00000000..e901d92e --- /dev/null +++ b/src/Rips_complex/doc/rips_complex_representation.png 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 @@ +<?xml version="1.0"?> +<!DOCTYPE ipe SYSTEM "ipe.dtd"> +<ipe version="70107" creator="Ipe 7.1.10"> +<info created="D:20150603143945" modified="D:20160928130224"/> +<ipestyle name="basic"> +<symbol name="arrow/arc(spx)"> +<path stroke="sym-stroke" fill="sym-stroke" pen="sym-pen"> +0 0 m +-1 0.333 l +-1 -0.333 l +h +</path> +</symbol> +<symbol name="arrow/farc(spx)"> +<path stroke="sym-stroke" fill="white" pen="sym-pen"> +0 0 m +-1 0.333 l +-1 -0.333 l +h +</path> +</symbol> +<symbol name="mark/circle(sx)" transformations="translations"> +<path fill="sym-stroke"> +0.6 0 0 0.6 0 0 e +0.4 0 0 0.4 0 0 e +</path> +</symbol> +<symbol name="mark/disk(sx)" transformations="translations"> +<path fill="sym-stroke"> +0.6 0 0 0.6 0 0 e +</path> +</symbol> +<symbol name="mark/fdisk(sfx)" transformations="translations"> +<group> +<path fill="sym-fill"> +0.5 0 0 0.5 0 0 e +</path> +<path fill="sym-stroke" fillrule="eofill"> +0.6 0 0 0.6 0 0 e +0.4 0 0 0.4 0 0 e +</path> +</group> +</symbol> +<symbol name="mark/box(sx)" transformations="translations"> +<path fill="sym-stroke" fillrule="eofill"> +-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 +</path> +</symbol> +<symbol name="mark/square(sx)" transformations="translations"> +<path fill="sym-stroke"> +-0.6 -0.6 m +0.6 -0.6 l +0.6 0.6 l +-0.6 0.6 l +h +</path> +</symbol> +<symbol name="mark/fsquare(sfx)" transformations="translations"> +<group> +<path fill="sym-fill"> +-0.5 -0.5 m +0.5 -0.5 l +0.5 0.5 l +-0.5 0.5 l +h +</path> +<path fill="sym-stroke" fillrule="eofill"> +-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 +</path> +</group> +</symbol> +<symbol name="mark/cross(sx)" transformations="translations"> +<group> +<path fill="sym-stroke"> +-0.43 -0.57 m +0.57 0.43 l +0.43 0.57 l +-0.57 -0.43 l +h +</path> +<path fill="sym-stroke"> +-0.43 0.57 m +0.57 -0.43 l +0.43 -0.57 l +-0.57 0.43 l +h +</path> +</group> +</symbol> +<symbol name="arrow/fnormal(spx)"> +<path stroke="sym-stroke" fill="white" pen="sym-pen"> +0 0 m +-1 0.333 l +-1 -0.333 l +h +</path> +</symbol> +<symbol name="arrow/pointed(spx)"> +<path stroke="sym-stroke" fill="sym-stroke" pen="sym-pen"> +0 0 m +-1 0.333 l +-0.8 0 l +-1 -0.333 l +h +</path> +</symbol> +<symbol name="arrow/fpointed(spx)"> +<path stroke="sym-stroke" fill="white" pen="sym-pen"> +0 0 m +-1 0.333 l +-0.8 0 l +-1 -0.333 l +h +</path> +</symbol> +<symbol name="arrow/linear(spx)"> +<path stroke="sym-stroke" pen="sym-pen"> +-1 0.333 m +0 0 l +-1 -0.333 l +</path> +</symbol> +<symbol name="arrow/fdouble(spx)"> +<path stroke="sym-stroke" fill="white" pen="sym-pen"> +0 0 m +-1 0.333 l +-1 -0.333 l +h +-1 0 m +-2 0.333 l +-2 -0.333 l +h +</path> +</symbol> +<symbol name="arrow/double(spx)"> +<path stroke="sym-stroke" fill="sym-stroke" pen="sym-pen"> +0 0 m +-1 0.333 l +-1 -0.333 l +h +-1 0 m +-2 0.333 l +-2 -0.333 l +h +</path> +</symbol> +<pen name="heavier" value="0.8"/> +<pen name="fat" value="1.2"/> +<pen name="ultrafat" value="2"/> +<symbolsize name="large" value="5"/> +<symbolsize name="small" value="2"/> +<symbolsize name="tiny" value="1.1"/> +<arrowsize name="large" value="10"/> +<arrowsize name="small" value="5"/> +<arrowsize name="tiny" value="3"/> +<color name="red" value="1 0 0"/> +<color name="green" value="0 1 0"/> +<color name="blue" value="0 0 1"/> +<color name="yellow" value="1 1 0"/> +<color name="orange" value="1 0.647 0"/> +<color name="gold" value="1 0.843 0"/> +<color name="purple" value="0.627 0.125 0.941"/> +<color name="gray" value="0.745"/> +<color name="brown" value="0.647 0.165 0.165"/> +<color name="navy" value="0 0 0.502"/> +<color name="pink" value="1 0.753 0.796"/> +<color name="seagreen" value="0.18 0.545 0.341"/> +<color name="turquoise" value="0.251 0.878 0.816"/> +<color name="violet" value="0.933 0.51 0.933"/> +<color name="darkblue" value="0 0 0.545"/> +<color name="darkcyan" value="0 0.545 0.545"/> +<color name="darkgray" value="0.663"/> +<color name="darkgreen" value="0 0.392 0"/> +<color name="darkmagenta" value="0.545 0 0.545"/> +<color name="darkorange" value="1 0.549 0"/> +<color name="darkred" value="0.545 0 0"/> +<color name="lightblue" value="0.678 0.847 0.902"/> +<color name="lightcyan" value="0.878 1 1"/> +<color name="lightgray" value="0.827"/> +<color name="lightgreen" value="0.565 0.933 0.565"/> +<color name="lightyellow" value="1 1 0.878"/> +<dashstyle name="dashed" value="[4] 0"/> +<dashstyle name="dotted" value="[1 3] 0"/> +<dashstyle name="dash dotted" value="[4 2 1 2] 0"/> +<dashstyle name="dash dot dotted" value="[4 2 1 2 1 2] 0"/> +<textsize name="large" value="\large"/> +<textsize name="small" value="\small"/> +<textsize name="tiny" value="\tiny"/> +<textsize name="Large" value="\Large"/> +<textsize name="LARGE" value="\LARGE"/> +<textsize name="huge" value="\huge"/> +<textsize name="Huge" value="\Huge"/> +<textsize name="footnote" value="\footnotesize"/> +<textstyle name="center" begin="\begin{center}" end="\end{center}"/> +<textstyle name="itemize" begin="\begin{itemize}" end="\end{itemize}"/> +<textstyle name="item" begin="\begin{itemize}\item{}" end="\end{itemize}"/> +<gridsize name="4 pts" value="4"/> +<gridsize name="8 pts (~3 mm)" value="8"/> +<gridsize name="16 pts (~6 mm)" value="16"/> +<gridsize name="32 pts (~12 mm)" value="32"/> +<gridsize name="10 pts (~3.5 mm)" value="10"/> +<gridsize name="20 pts (~7 mm)" value="20"/> +<gridsize name="14 pts (~5 mm)" value="14"/> +<gridsize name="28 pts (~10 mm)" value="28"/> +<gridsize name="56 pts (~20 mm)" value="56"/> +<anglesize name="90 deg" value="90"/> +<anglesize name="60 deg" value="60"/> +<anglesize name="45 deg" value="45"/> +<anglesize name="30 deg" value="30"/> +<anglesize name="22.5 deg" value="22.5"/> +<tiling name="falling" angle="-60" step="4" width="1"/> +<tiling name="rising" angle="30" step="4" width="1"/> +</ipestyle> +<page> +<layer name="alpha"/> +<view layers="alpha" active="alpha"/> +<path layer="alpha" matrix="1 0 0 1 0 -8" stroke="0"> +109.771 601.912 m +159.595 601.797 l +140.058 541.915 l +h +</path> +<path matrix="1 0 0 1 0 -8" stroke="0"> +79.8776 552.169 m +109.756 601.699 l +139.812 542.209 l +h +</path> +<path matrix="1 0 0 1 0.665417 -8.66542" stroke="0"> +69.8453 682.419 m +159.925 712.208 l +90.12 732.039 l +h +</path> +<text matrix="1 0 0 1 -230.178 14.1775" transformations="translations" pos="380 530" stroke="seagreen" type="label" width="98.916" height="8.307" depth="2.32" valign="baseline" size="large">One skeleton graph</text> +<text matrix="1 0 0 1 -212.333 10.6762" transformations="translations" pos="282.952 524.893" stroke="black" type="label" width="4.981" height="6.42" depth="0" valign="baseline">0</text> +<text matrix="1 0 0 1 -210.178 14.1775" transformations="translations" pos="352.708 510.349" stroke="black" type="label" width="4.981" height="6.42" depth="0" valign="baseline">1</text> +<text matrix="1 0 0 1 -210.178 14.1775" transformations="translations" pos="310.693 578.759" stroke="black" type="label" width="4.981" height="6.42" depth="0" valign="baseline">2</text> +<text matrix="1 0 0 1 -210.178 14.1775" transformations="translations" pos="375.332 578.49" stroke="black" type="label" width="4.981" height="6.42" depth="0" valign="baseline">3</text> +<text matrix="1 0 0 1 -210.178 14.1775" transformations="translations" pos="272.179 660.635" stroke="black" type="label" width="4.981" height="6.42" depth="0" valign="baseline">4</text> +<text matrix="1 0 0 1 -209.478 4.0238" transformations="translations" pos="296.419 724.197" stroke="black" type="label" width="4.981" height="6.42" depth="0" valign="baseline">5</text> +<text matrix="1 0 0 1 -210.178 14.1775" transformations="translations" pos="375.332 689.453" stroke="black" type="label" width="4.981" height="6.42" depth="0" valign="baseline">6</text> +<path matrix="1 0 0 1 30.6497 14.0396" stroke="black" pen="heavier"> +60 710 m +40 660 l +</path> +<path matrix="1 0 0 1 30.3739 13.9018" stroke="black" pen="heavier"> +40 660 m +130 690 l +</path> +<path matrix="1 0 0 1 29.8225 14.1775" stroke="black" pen="heavier"> +130 690 m +60 710 l +</path> +<path matrix="1 0 0 1 29.8225 14.1775" stroke="black" pen="heavier"> +40 660 m +80 580 l +</path> +<path matrix="1 0 0 1 29.8225 14.1775" stroke="black" pen="heavier"> +80 580 m +130 580 l +130 580 l +</path> +<path matrix="1 0 0 1 29.8225 14.1775" stroke="black" pen="heavier"> +130 580 m +110 520 l +</path> +<path matrix="1 0 0 1 29.8225 14.1775" stroke="black" pen="heavier"> +110 520 m +50 530 l +50 530 l +50 530 l +</path> +<path matrix="1 0 0 1 29.8225 14.1775" stroke="black" pen="heavier"> +50 530 m +80 580 l +</path> +<path matrix="1 0 0 1 29.8225 14.1775" stroke="black" pen="heavier"> +130 580 m +130 690 l +</path> +<use matrix="1 0 0 1 -209.478 4.0238" name="mark/fdisk(sfx)" pos="300 720" size="normal" stroke="black" fill="white"/> +<use matrix="1 0 0 1 -210.178 14.1775" name="mark/fdisk(sfx)" pos="280 660" size="normal" stroke="black" fill="white"/> +<use matrix="1 0 0 1 -210.178 14.1775" name="mark/fdisk(sfx)" pos="370 690" size="normal" stroke="black" fill="white"/> +<use matrix="1 0 0 1 -210.178 14.1775" name="mark/fdisk(sfx)" pos="370 580" size="normal" stroke="black" fill="white"/> +<use matrix="1 0 0 1 -210.178 14.1775" name="mark/fdisk(sfx)" pos="290 530" size="normal" stroke="black" fill="white"/> +<path matrix="1 0 0 1 -40 -16" stroke="black" pen="heavier"> +150.038 609.9 m +179.929 549.727 l +</path> +<use matrix="1 0 0 1 -210.178 14.1775" name="mark/fdisk(sfx)" pos="320 580" size="normal" stroke="black" fill="white"/> +<use matrix="1 0 0 1 -210.178 14.1775" name="mark/fdisk(sfx)" pos="350 520" size="normal" stroke="black" fill="white"/> +<path stroke="black" pen="heavier"> +158.7 593.269 m +81.4925 544.805 l +</path> +<path matrix="1 0 0 1 -17.9662 -17.9662" stroke="gray"> +256.324 639.958 m +370.055 639.958 l +</path> +<path matrix="1 0 0 1 -17.9662 -17.9662" stroke="gray"> +56.8567 0 0 56.8567 313.217 639.756 e +</path> +<use matrix="1 0 0 1 52.1387 -98.0941" name="mark/fdisk(sfx)" pos="300 720" size="normal" stroke="gray" fill="white"/> +<use matrix="1 0 0 1 -61.4926 -98.0942" name="mark/fdisk(sfx)" pos="300 720" size="normal" stroke="gray" fill="white"/> +<text matrix="1 0 0 1 -26.6167 -33.2708" transformations="translations" pos="295.735 657.944" stroke="gray" type="label" width="63.374" height="6.926" depth="1.93" valign="baseline">Rips threshold</text> +</page> +</ipe> diff --git a/src/Rips_complex/doc/rips_one_skeleton.png b/src/Rips_complex/doc/rips_one_skeleton.png Binary files differnew file mode 100644 index 00000000..1028770e --- /dev/null +++ b/src/Rips_complex/doc/rips_one_skeleton.png diff --git a/src/Rips_complex/example/CMakeLists.txt b/src/Rips_complex/example/CMakeLists.txt new file mode 100644 index 00000000..3962b15e --- /dev/null +++ b/src/Rips_complex/example/CMakeLists.txt @@ -0,0 +1,27 @@ +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}) + +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_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..2e63d9a6 --- /dev/null +++ b/src/Rips_complex/example/example_one_skeleton_rips_from_points.cpp @@ -0,0 +1,65 @@ +#include <gudhi/Rips_complex.h> +// to construct Rips_complex from a OFF file of points +#include <gudhi/Points_off_io.h> +// to construct a simplex_tree from rips complex +#include <gudhi/Simplex_tree.h> +#include <gudhi/distance_functions.h> + +#include <iostream> +#include <string> +#include <limits> // 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<double>; + 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>; + + std::vector<Point> 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_points(points, threshold, euclidean_distance<Filtration_value, Point>); + + Simplex_tree simplex; + if (rips_complex_from_points.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_off_file.cpp b/src/Rips_complex/example/example_rips_complex_from_off_file.cpp new file mode 100644 index 00000000..60050cea --- /dev/null +++ b/src/Rips_complex/example/example_rips_complex_from_off_file.cpp @@ -0,0 +1,72 @@ +#include <gudhi/Rips_complex.h> +// to construct Rips_complex from a OFF file of points +#include <gudhi/Points_off_io.h> +// to construct a simplex_tree from rips complex +#include <gudhi/Simplex_tree.h> +#include <gudhi/distance_functions.h> + +#include <iostream> +#include <string> + +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<float>; + using Simplex_tree = Gudhi::Simplex_tree<>; + using Filtration_value = Simplex_tree::Filtration_value; + using Rips_complex = Gudhi::rips_complex::Rips_complex<Filtration_value>; + + // ---------------------------------------------------------------------------- + // Init of a rips complex from an OFF file + // ---------------------------------------------------------------------------- + Gudhi::Points_off_reader<Point> off_reader(off_file_name); + Rips_complex rips_complex_from_file(off_reader.get_point_cloud(), threshold, + euclidean_distance<Filtration_value, Point>); + + 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/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/include/gudhi/Rips_complex.h b/src/Rips_complex/include/gudhi/Rips_complex.h new file mode 100644 index 00000000..aee6c969 --- /dev/null +++ b/src/Rips_complex/include/gudhi/Rips_complex.h @@ -0,0 +1,152 @@ +/* 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 <http://www.gnu.org/licenses/>. + */ + +#ifndef RIPS_COMPLEX_H_ +#define RIPS_COMPLEX_H_ + +#include <gudhi/Debug_utils.h> +#include <gudhi/graph_simplicial_complex.h> + +#include <boost/graph/adjacency_list.hpp> + +#include <iostream> +#include <vector> +#include <map> +#include <string> +#include <limits> // for numeric_limits +#include <utility> // for pair<> + + +namespace Gudhi { + +namespace rips_complex { + +/** + * \class Rips_complex + * \brief Rips complex data structure. + * + * \ingroup rips_complex + * + * \details + * 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. + * + * \tparam Filtration_value must meet `SimplicialComplexForRips` concept. + */ +template<typename Filtration_value> +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<typename InputPointRange, typename Point_d > + 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<Graph_t>::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 <typename SimplicialComplexForRips> + 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/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..f8be9748 --- /dev/null +++ b/src/Rips_complex/test/test_rips_complex.cpp @@ -0,0 +1,233 @@ +/* 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 <http://www.gnu.org/licenses/>. + */ + +#define BOOST_TEST_DYN_LINK +#define BOOST_TEST_MODULE "rips_complex" +#include <boost/test/unit_test.hpp> + +#include <cmath> // float comparison +#include <limits> +#include <string> +#include <vector> +#include <algorithm> // std::max + +#include <gudhi/Rips_complex.h> +// to construct Rips_complex from a OFF file of points +#include <gudhi/Points_off_io.h> +// to construct a simplex_tree from rips complex +#include <gudhi/Simplex_tree.h> +#include <gudhi/distance_functions.h> + +// Type definitions +using Point = std::vector<double>; +using Simplex_tree = Gudhi::Simplex_tree<>; +using Filtration_value = Simplex_tree::Filtration_value; +using Rips_complex = Gudhi::rips_complex::Rips_complex<Simplex_tree::Filtration_value>; + +bool are_almost_the_same(float a, float b) { + return std::fabs(a - b) < std::numeric_limits<float>::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<Point> off_reader(off_file_name); + Rips_complex rips_complex_from_file(off_reader.get_point_cloud(), rips_threshold, + euclidean_distance<Filtration_value, Point>); + + 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<Point> 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<double>(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<double>(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<Point>; + +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<double> 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<Point>); + + 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/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 <utility> // for pair #include <vector> -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<Graph_t, Edge_t, Filtration_value, Vertex_handle>(filegraph); // insert the graph in the simplex tree as 1-skeleton st.insert_graph(g); end = clock(); 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 <iostream> +#include <fstream> #include <string> #include <algorithm> #include <utility> // std::pair, std::make_pair @@ -19,19 +20,19 @@ using namespace Gudhi; typedef boost::mpl::list<Simplex_tree<>, Simplex_tree<Simplex_tree_options_fast_persistence>> 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<class typeST> 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<typename typeST::Simplex_handle, bool> typePairSimplexBool; - typedef std::vector<Vertex_handle> typeVectorVertex; + typedef std::vector<typename typeST::Vertex_handle> 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<class typeST, class typePairSimplexBool> 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<class typeST, class Filtration_value> @@ -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<typename typeST::Simplex_handle, bool> typePairSimplexBool; - typedef std::vector<Vertex_handle> typeVectorVertex; + typedef std::vector<typename typeST::Vertex_handle> typeVectorVertex; typedef std::pair<typeVectorVertex, Filtration_value> 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<typename typeST::Simplex_handle, bool> typePairSimplexBool; - typedef std::vector<Vertex_handle> typeVectorVertex; - typedef std::pair<typeVectorVertex, Filtration_value> typeSimplex; + typedef std::vector<typename typeST::Vertex_handle> typeVectorVertex; + typedef std::pair<typeVectorVertex, typename typeST::Filtration_value> 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<Vertex_handle>()); + std::sort(SimplexVector1.begin(), SimplexVector1.end(), std::greater<typename typeST::Vertex_handle>()); 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<Vertex_handle>()); + std::sort(SimplexVector2.begin(), SimplexVector2.end(), std::greater<typename typeST::Vertex_handle>()); 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<Vertex_handle>()); + std::sort(SimplexVector3.begin(), SimplexVector3.end(), std::greater<typename typeST::Vertex_handle>()); 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<Vertex_handle>()); + std::sort(SimplexVector5.begin(), SimplexVector5.end(), std::greater<typename typeST::Vertex_handle>()); 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<Vertex_handle>()); + std::sort(SimplexVector6.begin(), SimplexVector6.end(), std::greater<typename typeST::Vertex_handle>()); 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<Vertex_handle>& expected, int di } BOOST_AUTO_TEST_CASE_TEMPLATE(coface_on_simplex_tree, typeST, list_of_tested_variants) { - typedef std::vector<Vertex_handle> typeVectorVertex; + typedef std::vector<typename typeST::Vertex_handle> 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<Vertex_handle> simplex_result; + std::vector<typename typeST::Vertex_handle> simplex_result; std::vector<typename typeST::Simplex_handle> 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_handle> vertex = {3}; + std::vector<typename typeST::Vertex_handle> 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<Vertex_handle> typeVectorVertex; + typedef std::vector<typename typeST::Vertex_handle> 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 bb641b3c..59dd28e0 100644 --- a/src/Witness_complex/example/witness_complex_from_file.cpp +++ b/src/Witness_complex/example/witness_complex_from_file.cpp @@ -36,8 +36,9 @@ #include <string> #include <vector> -typedef std::vector< Vertex_handle > typeVectorVertex; +typedef std::vector< int > typeVectorVertex; typedef std::vector< std::vector <double> > Point_Vector; +typedef Gudhi::Simplex_tree<> Simplex_tree; int main(int argc, char * const argv[]) { if (argc != 3) { @@ -51,7 +52,7 @@ int main(int argc, char * const argv[]) { clock_t start, end; // Construct the Simplex Tree - Gudhi::Simplex_tree<> simplex_tree; + Simplex_tree simplex_tree; // Read the OFF file (input file name given as parameter) and triangulate points Gudhi::Points_off_reader<std::vector <double>> off_reader(off_file_name); @@ -69,7 +70,7 @@ int main(int argc, char * const argv[]) { std::vector<std::vector< int > > knn; Point_Vector landmarks; Gudhi::subsampling::pick_n_random_points(point_vector, 100, std::back_inserter(landmarks)); - Gudhi::witness_complex::construct_closest_landmark_table(point_vector, landmarks, knn); + Gudhi::witness_complex::construct_closest_landmark_table<Simplex_tree::Filtration_value>(point_vector, landmarks, knn); end = clock(); std::cout << "Landmark choice for " << nbL << " landmarks took " << static_cast<double>(end - start) / CLOCKS_PER_SEC << " s. \n"; diff --git a/src/Witness_complex/example/witness_complex_sphere.cpp b/src/Witness_complex/example/witness_complex_sphere.cpp index e6f88274..7ab86cc0 100644 --- a/src/Witness_complex/example/witness_complex_sphere.cpp +++ b/src/Witness_complex/example/witness_complex_sphere.cpp @@ -40,6 +40,8 @@ #include "generators.h" +typedef Gudhi::Simplex_tree<> Simplex_tree; + /** Write a gnuplot readable file. * Data range is a random access range of pairs (arg, value) */ @@ -62,7 +64,7 @@ int main(int argc, char * const argv[]) { clock_t start, end; // Construct the Simplex Tree - Gudhi::Simplex_tree<> simplex_tree; + Simplex_tree simplex_tree; std::vector< std::pair<int, double> > l_time; @@ -77,7 +79,7 @@ int main(int argc, char * const argv[]) { start = clock(); std::vector<std::vector< int > > knn; Gudhi::subsampling::pick_n_random_points(point_vector, 100, std::back_inserter(landmarks)); - Gudhi::witness_complex::construct_closest_landmark_table(point_vector, landmarks, knn); + Gudhi::witness_complex::construct_closest_landmark_table<Simplex_tree::Filtration_value>(point_vector, landmarks, knn); // Compute witness complex Gudhi::witness_complex::witness_complex(knn, number_of_landmarks, point_vector[0].size(), simplex_tree); diff --git a/src/Witness_complex/include/gudhi/Construct_closest_landmark_table.h b/src/Witness_complex/include/gudhi/Construct_closest_landmark_table.h index ef711c34..ec93ae71 100644 --- a/src/Witness_complex/include/gudhi/Construct_closest_landmark_table.h +++ b/src/Witness_complex/include/gudhi/Construct_closest_landmark_table.h @@ -51,7 +51,8 @@ namespace witness_complex { * Closest_landmark_range needs to have push_back operation. */ - template <typename WitnessContainer, + template <typename FiltrationValue, + typename WitnessContainer, typename LandmarkContainer, typename KNearestNeighbours> void construct_closest_landmark_table(WitnessContainer const &points, @@ -72,7 +73,8 @@ namespace witness_complex { int landmarks_i = 0; for (landmarks_it = landmarks.begin(), landmarks_i = 0; landmarks_it != landmarks.end(); ++landmarks_it, landmarks_i++) { - dist_i dist = std::make_pair(euclidean_distance(points[points_i], *landmarks_it), landmarks_i); + dist_i dist = std::make_pair(euclidean_distance<FiltrationValue>(points[points_i], *landmarks_it), + landmarks_i); l_heap.push(dist); } for (int i = 0; i < dim + 1; i++) { diff --git a/src/Witness_complex/include/gudhi/Landmark_choice_by_furthest_point.h b/src/Witness_complex/include/gudhi/Landmark_choice_by_furthest_point.h new file mode 100644 index 00000000..bcb89e00 --- /dev/null +++ b/src/Witness_complex/include/gudhi/Landmark_choice_by_furthest_point.h @@ -0,0 +1,105 @@ +/* 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): Siargey Kachanovich + * + * Copyright (C) 2015 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/>. + */ + +#ifndef LANDMARK_CHOICE_BY_FURTHEST_POINT_H_ +#define LANDMARK_CHOICE_BY_FURTHEST_POINT_H_ + +#include <boost/range/size.hpp> + +#include <limits> // for numeric_limits<> +#include <iterator> +#include <algorithm> // for sort +#include <vector> + +namespace Gudhi { + +namespace witness_complex { + + typedef std::vector<int> typeVectorVertex; + + /** + * \ingroup witness_complex + * \brief Landmark choice strategy by iteratively adding the furthest witness from the + * current landmark set as the new landmark. + * \details It chooses nbL landmarks from a random access range `points` and + * writes {witness}*{closest landmarks} matrix in `knn`. + * + * The type KNearestNeighbors can be seen as + * Witness_range<Closest_landmark_range<Vertex_handle>>, where + * Witness_range and Closest_landmark_range are random access ranges + * + */ + + template <typename KNearestNeighbours, + typename Point_random_access_range> + void landmark_choice_by_furthest_point(Point_random_access_range const &points, + int nbL, + KNearestNeighbours &knn) { + int nb_points = boost::size(points); + assert(nb_points >= nbL); + // distance matrix witness x landmarks + std::vector<std::vector<double>> wit_land_dist(nb_points, std::vector<double>()); + // landmark list + typeVectorVertex chosen_landmarks; + + knn = KNearestNeighbours(nb_points, std::vector<int>()); + int current_number_of_landmarks = 0; // counter for landmarks + double curr_max_dist = 0; // used for defining the furhest point from L + const double infty = std::numeric_limits<double>::infinity(); // infinity (see next entry) + std::vector< double > dist_to_L(nb_points, infty); // vector of current distances to L from points + + // TODO(SK) Consider using rand_r(...) instead of rand(...) for improved thread safety + // or better yet std::uniform_int_distribution + int rand_int = rand() % nb_points; + int curr_max_w = rand_int; // For testing purposes a pseudo-random number is used here + + for (current_number_of_landmarks = 0; current_number_of_landmarks != nbL; current_number_of_landmarks++) { + // curr_max_w at this point is the next landmark + chosen_landmarks.push_back(curr_max_w); + unsigned i = 0; + for (auto& p : points) { + double curr_dist = euclidean_distance<double>(p, *(std::begin(points) + chosen_landmarks[current_number_of_landmarks])); + wit_land_dist[i].push_back(curr_dist); + knn[i].push_back(current_number_of_landmarks); + if (curr_dist < dist_to_L[i]) + dist_to_L[i] = curr_dist; + ++i; + } + curr_max_dist = 0; + for (i = 0; i < dist_to_L.size(); i++) + if (dist_to_L[i] > curr_max_dist) { + curr_max_dist = dist_to_L[i]; + curr_max_w = i; + } + } + for (int i = 0; i < nb_points; ++i) + std::sort(std::begin(knn[i]), + std::end(knn[i]), + [&wit_land_dist, i](int a, int b) { + return wit_land_dist[i][a] < wit_land_dist[i][b]; }); + } + +} // namespace witness_complex + +} // namespace Gudhi + +#endif // LANDMARK_CHOICE_BY_FURTHEST_POINT_H_ diff --git a/src/Witness_complex/include/gudhi/Landmark_choice_by_random_point.h b/src/Witness_complex/include/gudhi/Landmark_choice_by_random_point.h new file mode 100644 index 00000000..b5aab9d5 --- /dev/null +++ b/src/Witness_complex/include/gudhi/Landmark_choice_by_random_point.h @@ -0,0 +1,96 @@ +/* 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): Siargey Kachanovich + * + * Copyright (C) 2015 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/>. + */ + +#ifndef LANDMARK_CHOICE_BY_RANDOM_POINT_H_ +#define LANDMARK_CHOICE_BY_RANDOM_POINT_H_ + +#include <boost/range/size.hpp> + +#include <queue> // for priority_queue<> +#include <utility> // for pair<> +#include <iterator> +#include <vector> +#include <set> + +namespace Gudhi { + +namespace witness_complex { + + /** + * \ingroup witness_complex + * \brief Landmark choice strategy by taking random vertices for landmarks. + * \details It chooses nbL distinct landmarks from a random access range `points` + * and outputs a matrix {witness}*{closest landmarks} in knn. + * + * The type KNearestNeighbors can be seen as + * Witness_range<Closest_landmark_range<Vertex_handle>>, where + * Witness_range and Closest_landmark_range are random access ranges and + * Vertex_handle is the label type of a vertex in a simplicial complex. + * Closest_landmark_range needs to have push_back operation. + */ + + template <typename KNearestNeighbours, + typename Point_random_access_range> + void landmark_choice_by_random_point(Point_random_access_range const &points, + int nbL, + KNearestNeighbours &knn) { + int nbP = boost::size(points); + assert(nbP >= nbL); + std::set<int> landmarks; + int current_number_of_landmarks = 0; // counter for landmarks + + // TODO(SK) Consider using rand_r(...) instead of rand(...) for improved thread safety + int chosen_landmark = rand() % nbP; + for (current_number_of_landmarks = 0; current_number_of_landmarks != nbL; current_number_of_landmarks++) { + while (landmarks.find(chosen_landmark) != landmarks.end()) + chosen_landmark = rand() % nbP; + landmarks.insert(chosen_landmark); + } + + int dim = boost::size(*std::begin(points)); + typedef std::pair<double, int> dist_i; + typedef bool (*comp)(dist_i, dist_i); + knn = KNearestNeighbours(nbP); + for (int points_i = 0; points_i < nbP; points_i++) { + std::priority_queue<dist_i, std::vector<dist_i>, comp> l_heap([](dist_i j1, dist_i j2) { + return j1.first > j2.first; + }); + std::set<int>::iterator landmarks_it; + int landmarks_i = 0; + for (landmarks_it = landmarks.begin(), landmarks_i = 0; landmarks_it != landmarks.end(); + ++landmarks_it, landmarks_i++) { + dist_i dist = std::make_pair(euclidean_distance<double>(points[points_i], points[*landmarks_it]), landmarks_i); + l_heap.push(dist); + } + for (int i = 0; i < dim + 1; i++) { + dist_i dist = l_heap.top(); + knn[points_i].push_back(dist.second); + l_heap.pop(); + } + } + } + +} // namespace witness_complex + +} // namespace Gudhi + +#endif // LANDMARK_CHOICE_BY_RANDOM_POINT_H_ 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 <vector> typedef Gudhi::Simplex_tree<> Simplex_tree; -typedef std::vector< Vertex_handle > typeVectorVertex; +typedef std::vector< int > typeVectorVertex; typedef Gudhi::witness_complex::Witness_complex<Simplex_tree> 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 d40bbf14..b7067f87 100644 --- a/src/Witness_complex/test/witness_complex_points.cpp +++ b/src/Witness_complex/test/witness_complex_points.cpp @@ -34,7 +34,7 @@ #include <vector> typedef std::vector<double> Point; -typedef std::vector< Vertex_handle > typeVectorVertex; +typedef std::vector< int > typeVectorVertex; typedef Gudhi::Simplex_tree<> Simplex_tree; typedef Gudhi::witness_complex::Witness_complex<Simplex_tree> WitnessComplex; @@ -51,7 +51,7 @@ BOOST_AUTO_TEST_CASE(witness_complex_points) { // First test: random choice Simplex_tree complex1; Gudhi::subsampling::pick_n_random_points(points, 100, std::back_inserter(landmarks)); - Gudhi::witness_complex::construct_closest_landmark_table(points, landmarks, knn); + Gudhi::witness_complex::construct_closest_landmark_table<Simplex_tree::Filtration_value>(points, landmarks, knn); assert(!knn.empty()); WitnessComplex witnessComplex1(knn, 100, 3, complex1); BOOST_CHECK(witnessComplex1.is_witness_complex(knn, b_print_output)); diff --git a/src/cmake/modules/GUDHI_user_version_target.txt b/src/cmake/modules/GUDHI_user_version_target.txt index 51553e7e..7542629e 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;Spatial_searching;Subsampling;Tangential_complex;Witness_complex") + set(GUDHI_MODULES "common;Alpha_complex;Bitmap_cubical_complex;Contraction;Hasse_complex;Persistent_cohomology;Rips_complex;Simplex_tree;Skeleton_blocker;Spatial_searching;Subsampling;Tangential_complex;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 1a2cb6ba..85310f66 100644 --- a/src/common/doc/main_page.h +++ b/src/common/doc/main_page.h @@ -55,6 +55,23 @@ <b>User manual:</b> \ref cubical_complex - <b>Reference manual:</b> Gudhi::cubical_complex::Bitmap_cubical_complex </td> </tr> + \subsection RipsComplexDataStructure Rips complex + \image html "rips_complex_representation.png" "Rips complex representation" +<table border="0"> + <tr> + <td width="25%"> + <b>Author:</b> Clément Maria<br> + <b>Introduced in:</b> GUDHI 1.4.0<br> + <b>Copyright:</b> GPL v3<br> + </td> + <td width="75%"> + Rips_complex is a simplicial complex constructed from a one skeleton graph.<br> + The filtration value of each edge is computed from a user-given distance function and is inserted until a + user-given threshold value.<br> + <b>User manual:</b> \ref rips_complex - <b>Reference manual:</b> Gudhi::rips_complex::Rips_complex + </td> + </tr> +</table> </table> \subsection SimplexTreeDataStructure Simplex tree \image html "Simplex_tree_representation.png" "Simplex tree representation" @@ -202,8 +219,6 @@ make \endverbatim * * The following example requires the <a target="_blank" href="http://gmplib.org/">GNU Multiple Precision Arithmetic * Library</a> (GMP) and will not be built if GMP is not installed: - * \li <a href="_persistent_cohomology_2performance_rips_persistence_8cpp-example.html"> - * Persistent_cohomology/performance_rips_persistence.cpp</a> * \li <a href="_persistent_cohomology_2rips_multifield_persistence_8cpp-example.html"> * Persistent_cohomology/rips_multifield_persistence.cpp</a> * @@ -291,8 +306,6 @@ make \endverbatim * Persistent_cohomology/alpha_complex_persistence.cpp</a> * \li <a href="_persistent_cohomology_2rips_persistence_via_boundary_matrix_8cpp-example.html"> * Persistent_cohomology/rips_persistence_via_boundary_matrix.cpp</a> - * \li <a href="_persistent_cohomology_2performance_rips_persistence_8cpp-example.html"> - * Persistent_cohomology/performance_rips_persistence.cpp</a> * \li <a href="_persistent_cohomology_2persistence_from_file_8cpp-example.html"> * Persistent_cohomology/persistence_from_file.cpp</a> * \li <a href="_persistent_cohomology_2persistence_from_simple_simplex_tree_8cpp-example.html"> @@ -307,6 +320,10 @@ make \endverbatim * Persistent_cohomology/periodic_alpha_complex_3d_persistence.cpp</a> * \li <a href="_persistent_cohomology_2custom_persistence_sort_8cpp-example.html"> * Persistent_cohomology/custom_persistence_sort.cpp</a> + * \li <a href="_rips_complex_2example_one_skeleton_rips_from_points_8cpp-example.html"> + * Rips_complex/example_one_skeleton_rips_from_points.cpp</a> + * \li <a href="_rips_complex_2example_rips_complex_from_off_file_8cpp-example.html"> + * Rips_complex/example_rips_complex_from_off_file.cpp</a> * * \section Contributions Bug reports and contributions * Please help us improving the quality of the GUDHI library. You may report bugs or suggestions to: @@ -341,7 +358,6 @@ make \endverbatim * @example Persistent_cohomology/alpha_complex_3d_persistence.cpp * @example Persistent_cohomology/alpha_complex_persistence.cpp * @example Persistent_cohomology/rips_persistence_via_boundary_matrix.cpp - * @example Persistent_cohomology/performance_rips_persistence.cpp * @example Persistent_cohomology/periodic_alpha_complex_3d_persistence.cpp * @example Persistent_cohomology/persistence_from_file.cpp * @example Persistent_cohomology/persistence_from_simple_simplex_tree.cpp @@ -349,6 +365,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/example_alpha_shapes_3_simplex_tree_from_off_file.cpp diff --git a/src/common/include/gudhi/distance_functions.h b/src/common/include/gudhi/distance_functions.h index cd518581..006df275 100644 --- a/src/common/include/gudhi/distance_functions.h +++ b/src/common/include/gudhi/distance_functions.h @@ -28,13 +28,13 @@ /* Compute the Euclidean distance between two Points given * 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 dist = 0.; +template< typename Filtration_value, typename Point > +Filtration_value euclidean_distance(const Point &p1,const Point &p2) { + Filtration_value dist = 0.; auto it1 = p1.begin(); auto it2 = p2.begin(); for (; it1 != p1.end(); ++it1, ++it2) { - double tmp = *it1 - *it2; + Filtration_value tmp = (double)(*it1) - (double)(*it2); dist += tmp*tmp; } return std::sqrt(dist); 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_ 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<Graph_t>::vertex_iterator vi, vi_end; + typename boost::graph_traits<Graph_t>::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); |
