From 8837fea64910e8f2e45bebe25c3733a8250ebca8 Mon Sep 17 00:00:00 2001 From: skachano Date: Thu, 18 Feb 2016 11:13:57 +0000 Subject: Improved the benching program a bit git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/relaxed-witness@1030 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: dc3d13fa785bee42a574134519b05e98b934d6f5 --- src/Persistent_cohomology/example/CMakeLists.txt | 16 ++ .../example/relaxed_witness_persistence.cpp | 233 ++++++++++++++++++--- .../include/gudhi/Relaxed_witness_complex.h | 18 +- 3 files changed, 229 insertions(+), 38 deletions(-) diff --git a/src/Persistent_cohomology/example/CMakeLists.txt b/src/Persistent_cohomology/example/CMakeLists.txt index 95506631..772193c4 100644 --- a/src/Persistent_cohomology/example/CMakeLists.txt +++ b/src/Persistent_cohomology/example/CMakeLists.txt @@ -31,6 +31,7 @@ add_test(persistence_from_file_3_3_100 ${CMAKE_CURRENT_BINARY_DIR}/persistence_f if(GMPXX_FOUND AND GMP_FOUND) message("GMPXX_LIBRARIES = ${GMPXX_LIBRARIES}") message("GMP_LIBRARIES = ${GMP_LIBRARIES}") + message("CGAL_LIBRARY = ${CGAL_LIBRARY}") 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}) @@ -40,6 +41,21 @@ if(GMPXX_FOUND AND GMP_FOUND) target_link_libraries(performance_rips_persistence ${Boost_SYSTEM_LIBRARY} ${Boost_PROGRAM_OPTIONS_LIBRARY} ${GMPXX_LIBRARIES} ${GMP_LIBRARIES}) if(CGAL_FOUND) + include( ${CGAL_USE_FILE} ) + # In CMakeLists.txt, when include(${CGAL_USE_FILE}), CXX_FLAGS are overwritten. + # cf. http://doc.cgal.org/latest/Manual/installation.html#title40 + # A workaround is to add "-std=c++11" again. + # A fix would be to use https://cmake.org/cmake/help/v3.1/prop_gbl/CMAKE_CXX_KNOWN_FEATURES.html + # or even better https://cmake.org/cmake/help/v3.1/variable/CMAKE_CXX_STANDARD.html + # but it implies to use cmake version 3.1 at least. + if(NOT MSVC) + include(CheckCXXCompilerFlag) + CHECK_CXX_COMPILER_FLAG(-std=c++11 COMPILER_SUPPORTS_CXX11) + if(COMPILER_SUPPORTS_CXX11) + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++11") + endif() + endif() + # - End of workaround if (CMAKE_BUILD_TYPE MATCHES Debug) # For programs to be more verbose add_definitions(-DDEBUG_TRACES) diff --git a/src/Witness_complex/example/relaxed_witness_persistence.cpp b/src/Witness_complex/example/relaxed_witness_persistence.cpp index 09fe0d31..b8e4866f 100644 --- a/src/Witness_complex/example/relaxed_witness_persistence.cpp +++ b/src/Witness_complex/example/relaxed_witness_persistence.cpp @@ -50,41 +50,63 @@ using namespace Gudhi::persistent_cohomology; typedef std::vector Point_Vector; typedef Relaxed_witness_complex< Simplex_tree<> > RelaxedWitnessComplex; -int main (int argc, char * const argv[]) -{ - if (argc != 6) { - std::cerr << "Usage: " << argv[0] - << " nbP nbL dim alpha limD\n"; - return 0; +/** + * \brief Customized version of read_points + * which takes into account a possible nbP first line + * + */ +inline void +read_points_cust(std::string file_name, std::vector< std::vector< double > > & points) { + std::ifstream in_file(file_name.c_str(), std::ios::in); + if (!in_file.is_open()) { + std::cerr << "Unable to open file " << file_name << std::endl; + return; + } + std::string line; + double x; + while (getline(in_file, line)) { + std::vector< double > point; + std::istringstream iss(line); + while (iss >> x) { + point.push_back(x); } - /* - boost::filesystem::path p; - for (; argc > 2; --argc, ++argv) - p /= argv[1]; - */ - - int nbP = atoi(argv[1]); - int nbL = atoi(argv[2]); - int dim = atoi(argv[3]); - double alpha = atof(argv[4]); - int limD = atoi(argv[5]); - //Construct the Simplex Tree - clock_t start, end; + if (point.size() != 1) + points.push_back(point); + } + in_file.close(); +} - // Construct the Simplex Tree - Simplex_tree<> simplex_tree; +void output_experiment_information(char * const file_name) +{ + std::cout << "Enter a valid experiment number. Usage: " + << file_name << " exp_no options\n"; + std::cout << "Experiment description:\n" + << "0 nbP nbL dim alpha limD: " + << "Build persistence diagram on relaxed witness complex " + << "built from a point cloud on (dim-1)-dimensional sphere " + << "consisting of nbP witnesses and nbL landmarks. " + << "The maximal relaxation is alpha and the limit on simplicial complex " + << "dimension is limD.\n"; + std::cout << "1 file_name nbL alpha limD: " + << "Build persistence diagram on relaxed witness complex " + << "build from a point cloud stored in a file and nbL landmarks. " + << "The maximal relaxation is alpha and the limit on simplicial complex dimension is limD\n"; +} - // Read the point file - Point_Vector point_vector; - generate_points_sphere(point_vector, nbP, dim); - std::cout << "Successfully generated " << point_vector.size() << " points.\n"; - std::cout << "Ambient dimension is " << point_vector[0].size() << ".\n"; +void rw_experiment(Point_Vector & point_vector, int nbL, FT alpha, int limD) +{ + clock_t start, end; + Simplex_tree<> simplex_tree; // Choose landmarks std::vector > knn; std::vector > distances; + start = clock(); Gudhi::witness_complex::landmark_choice_by_random_knn(point_vector, nbL, alpha, limD, knn, distances); - + end = clock(); + double time = static_cast(end - start) / CLOCKS_PER_SEC; + std::cout << "Choice of " << nbL << " landmarks took " + << time << " s. \n"; // Compute witness complex start = clock(); RelaxedWitnessComplex(distances, @@ -94,15 +116,168 @@ int main (int argc, char * const argv[]) alpha, limD); end = clock(); - double time = static_cast(end - start) / CLOCKS_PER_SEC; + time = static_cast(end - start) / CLOCKS_PER_SEC; std::cout << "Witness complex for " << nbL << " landmarks took " << time << " s. \n"; + std::cout << "The complex contains " << simplex_tree.num_simplices() << " simplices \n"; // std::cout << simplex_tree << "\n"; // Compute the persistence diagram of the complex persistent_cohomology::Persistent_cohomology< Simplex_tree<>, Field_Zp > pcoh(simplex_tree, false); int p = 3; pcoh.init_coefficients( p ); //initilizes the coefficient field for homology - pcoh.compute_persistent_cohomology( alpha/10 ); + start = clock(); + pcoh.compute_persistent_cohomology( alpha/5 ); + end = clock(); + time = static_cast(end - start) / CLOCKS_PER_SEC; + std::cout << "Persistence diagram took " + << time << " s. \n"; pcoh.output_diagram(); } + +void rips_experiment(Point_Vector & points, double threshold, int dim_max) +{ + typedef std::vector Point_t; + typedef Simplex_tree ST; + clock_t start, end; + ST st; + + // Compute the proximity graph of the points + start = clock(); + Graph_t prox_graph = compute_proximity_graph(points, threshold + , euclidean_distance); + // Construct the Rips complex in a 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); + end = clock(); + + double time = static_cast(end - start) / CLOCKS_PER_SEC; + std::cout << "Rips complex took " + << time << " s. \n"; + std::cout << "The complex contains " << st.num_simplices() << " simplices \n"; + //std::cout << " and has dimension " << st.dimension() << " \n"; + + // Sort the simplices in the order of the filtration + st.initialize_filtration(); + + // Compute the persistence diagram of the complex + persistent_cohomology::Persistent_cohomology pcoh(st); + // initializes the coefficient field for homology + int p = 3; + double min_persistence = threshold/5; + pcoh.init_coefficients(p); + pcoh.compute_persistent_cohomology(min_persistence); + pcoh.output_diagram(); +} + +int experiment0 (int argc, char * const argv[]) +{ + if (argc != 7) { + std::cerr << "Usage: " << argv[0] + << " 0 nbP nbL dim alpha limD\n"; + return 0; + } + /* + boost::filesystem::path p; + for (; argc > 2; --argc, ++argv) + p /= argv[1]; + */ + + int nbP = atoi(argv[2]); + int nbL = atoi(argv[3]); + int dim = atoi(argv[4]); + double alpha = atof(argv[5]); + int limD = atoi(argv[6]); + + // Read the point file + Point_Vector point_vector; + generate_points_sphere(point_vector, nbP, dim); + std::cout << "Successfully generated " << point_vector.size() << " points.\n"; + std::cout << "Ambient dimension is " << point_vector[0].size() << ".\n"; + + rw_experiment(point_vector, nbL, alpha, limD); +} + +int experiment1 (int argc, char * const argv[]) +{ + if (argc != 3) { + std::cerr << "Usage: " << argv[0] + << " 1 file_name\n"; + return 0; + } + /* + boost::filesystem::path p; + for (; argc > 2; --argc, ++argv) + p /= argv[1]; + */ + + std::string file_name = argv[2]; + + // Read the point file + Point_Vector point_vector; + read_points_cust(file_name, point_vector); + std::cout << "The file contains " << point_vector.size() << " points.\n"; + std::cout << "Ambient dimension is " << point_vector[0].size() << ".\n"; + + bool ok = false; + int nbL, limD; + double alpha; + while (!ok) { + std::cout << "Relaxed witness complex: parameters nbL, alpha, limD.\n"; + std::cout << "Enter nbL: "; + std::cin >> nbL; + std::cout << "Enter alpha: "; + std::cin >> alpha; + std::cout << "Enter limD: "; + std::cin >> limD; + std::cout << "Start relaxed witness complex...\n"; + rw_experiment(point_vector, nbL, alpha, limD); + std::cout << "Is the result correct? [y/n]: "; + char answer; + std::cin >> answer; + switch (answer) { + case 'n': + ok = false; break; + default : + ok = true; break; + } + } + ok = false; + while (!ok) { + std::cout << "Rips complex: parameters threshold, limD.\n"; + std::cout << "Enter threshold: "; + std::cin >> alpha; + std::cout << "Enter limD: "; + std::cin >> limD; + std::cout << "Start Rips complex...\n"; + rips_experiment(point_vector, alpha, limD); + std::cout << "Is the result correct? [y/n]: "; + char answer; + std::cin >> answer; + switch (answer) { + case 'n': + ok = false; break; + default : + ok = true; break; + } + } +} + +int main (int argc, char * const argv[]) +{ + if (argc == 1) { + output_experiment_information(argv[0]); + return 1; + } + switch (atoi(argv[1])) { + case 0 : + return experiment0(argc, argv); + case 1 : + return experiment1(argc, argv); + default : + output_experiment_information(argv[0]); + return 1; + } +} diff --git a/src/Witness_complex/include/gudhi/Relaxed_witness_complex.h b/src/Witness_complex/include/gudhi/Relaxed_witness_complex.h index 6fc71a05..04dc37d4 100644 --- a/src/Witness_complex/include/gudhi/Relaxed_witness_complex.h +++ b/src/Witness_complex/include/gudhi/Relaxed_witness_complex.h @@ -70,17 +70,17 @@ namespace witness_complex { template< class Simplicial_complex > class Relaxed_witness_complex { - private: - struct Active_witness { - int witness_id; - int landmark_id; - - Active_witness(int witness_id_, int landmark_id_) - : witness_id(witness_id_), - landmark_id(landmark_id_) { } +private: + struct Active_witness { + int witness_id; + int landmark_id; + + Active_witness(int witness_id_, int landmark_id_) + : witness_id(witness_id_), + landmark_id(landmark_id_) { } }; - private: +private: typedef typename Simplicial_complex::Simplex_handle Simplex_handle; typedef typename Simplicial_complex::Vertex_handle Vertex_handle; -- cgit v1.2.3