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 --- .../example/relaxed_witness_persistence.cpp | 233 ++++++++++++++++++--- 1 file changed, 204 insertions(+), 29 deletions(-) (limited to 'src/Witness_complex/example') 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; + } +} -- cgit v1.2.3