From 1c09c91ddf4d38196a91bbff5ae432fb13be6762 Mon Sep 17 00:00:00 2001 From: skachano Date: Wed, 28 Sep 2016 11:54:31 +0000 Subject: All work up until now git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/relaxed-witness@1578 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: a9379060f1edfbc419e259b2e207ebc608798803 --- src/Witness_complex/example/CMakeLists.txt | 2 + .../example/Landmark_choice_random_knn.h | 2 +- .../example/Landmark_choice_sparsification.h | 6 +- src/Witness_complex/example/a0_persistence.cpp | 8 +- src/Witness_complex/example/bench_rwit.cpp | 329 +++++++++++------ src/Witness_complex/example/generators.h | 17 +- src/Witness_complex/example/output.h | 9 +- src/Witness_complex/include/gudhi/Good_links.h | 399 ++++++++++++++++++++- 8 files changed, 632 insertions(+), 140 deletions(-) (limited to 'src') diff --git a/src/Witness_complex/example/CMakeLists.txt b/src/Witness_complex/example/CMakeLists.txt index 80396392..d6a958b7 100644 --- a/src/Witness_complex/example/CMakeLists.txt +++ b/src/Witness_complex/example/CMakeLists.txt @@ -42,10 +42,12 @@ if(CGAL_FOUND) target_link_libraries(a0_persistence ${Boost_SYSTEM_LIBRARY} ${CGAL_LIBRARY}) add_test( a0_persistence_10 ${CMAKE_CURRENT_BINARY_DIR}/a0_persistence 10) add_executable ( bench_rwit bench_rwit.cpp ) + target_link_libraries(bench_rwit ${Boost_SYSTEM_LIBRARY} ${Boost_PROGRAM_OPTIONS_LIBRARY}) add_executable ( relaxed_delaunay relaxed_delaunay.cpp ) target_link_libraries(relaxed_delaunay ${Boost_SYSTEM_LIBRARY} ${CGAL_LIBRARY}) add_test( relaxed_delaunay_10 ${CMAKE_CURRENT_BINARY_DIR}/relaxed_delaunay 10) + add_executable ( generate_torus generate_torus.cpp ) else() message(WARNING "Eigen3 not found. Version 3.1.0 is required for witness_complex_sphere example.") endif() diff --git a/src/Witness_complex/example/Landmark_choice_random_knn.h b/src/Witness_complex/example/Landmark_choice_random_knn.h index 3797b4c5..fb91e116 100644 --- a/src/Witness_complex/example/Landmark_choice_random_knn.h +++ b/src/Witness_complex/example/Landmark_choice_random_knn.h @@ -131,7 +131,7 @@ typedef Neighbor_search::Tree Tree; distances[points_i].push_back(sqrt(search_it->second)); knn[points_i].push_back((search_it++)->first); } - std::cout << "k = " << knn[points_i].size() << std::endl; + //std::cout << "k = " << knn[points_i].size() << std::endl; } } diff --git a/src/Witness_complex/example/Landmark_choice_sparsification.h b/src/Witness_complex/example/Landmark_choice_sparsification.h index 1052b0c4..d5d5fba1 100644 --- a/src/Witness_complex/example/Landmark_choice_sparsification.h +++ b/src/Witness_complex/example/Landmark_choice_sparsification.h @@ -74,7 +74,7 @@ namespace witness_complex { FT mu_epsilon, Point_random_access_range & landmarks) { - int nbP = points.end() - points.begin(); + unsigned nbP = points.end() - points.begin(); assert(nbP >= nbL); CGAL::Random rand; // TODO(SK) Consider using rand_r(...) instead of rand(...) for improved thread safety @@ -87,7 +87,7 @@ namespace witness_complex { typename Tree::Splitter(), points_traits); - for (int points_i = 0; points_i < nbP; points_i++) { + for (unsigned points_i = 0; points_i < nbP; points_i++) { if (dropped_points[points_i]) continue; Point_d & w = points[points_i]; @@ -98,7 +98,7 @@ namespace witness_complex { dropped_points[i] = true; } - for (int points_i = 0; points_i < nbP; points_i++) { + for (unsigned points_i = 0; points_i < nbP; points_i++) { if (dropped_points[points_i]) landmarks.push_back(points[points_i]); } diff --git a/src/Witness_complex/example/a0_persistence.cpp b/src/Witness_complex/example/a0_persistence.cpp index 295c7115..422185ca 100644 --- a/src/Witness_complex/example/a0_persistence.cpp +++ b/src/Witness_complex/example/a0_persistence.cpp @@ -138,6 +138,12 @@ void rw_experiment(Point_Vector & point_vector, int nbL, FT alpha2, int limD, FT std::cout << "Witness complex for " << nbL << " landmarks took " << time << " s. \n"; std::cout << "The complex contains " << simplex_tree.num_simplices() << " simplices \n"; + int streedim = 0; + for (auto s: simplex_tree.complex_simplex_range()) + if (simplex_tree.dimension(s) > streedim) + streedim = simplex_tree.dimension(s); + std::cout << "Dimension of the simplicial complex is " << streedim << std::endl; + //std::cout << simplex_tree << "\n"; // Compute the persistence diagram of the complex @@ -348,7 +354,7 @@ int experiment1 (int argc, char * const argv[]) 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"; - + Simplex_tree<> simplex_tree; std::vector > knn; std::vector > distances; diff --git a/src/Witness_complex/example/bench_rwit.cpp b/src/Witness_complex/example/bench_rwit.cpp index f83f9f42..2d3a009c 100644 --- a/src/Witness_complex/example/bench_rwit.cpp +++ b/src/Witness_complex/example/bench_rwit.cpp @@ -26,6 +26,7 @@ #include #include #include +#include #include "Landmark_choice_random_knn.h" #include "Landmark_choice_sparsification.h" @@ -44,6 +45,8 @@ #include #include +#include + #include "generators.h" #include "output.h" #include "output_tikz.h" @@ -53,14 +56,96 @@ using namespace Gudhi::witness_complex; using namespace Gudhi::persistent_cohomology; typedef std::vector Point_Vector; -typedef Simplex_tree STree; -//typedef Simplex_tree<> STree; +//typedef Simplex_tree STree; +typedef Simplex_tree<> STree; typedef A0_complex< STree> A0Complex; typedef STree::Simplex_handle Simplex_handle; typedef A0_complex< STree > SRWit; typedef Relaxed_witness_complex< STree > WRWit; +/** Program options *************************************************************** +*********************************************************************************** +*********************************************************************************** +*********************************************************************************** +**********************************************************************************/ + + +void program_options(int argc, char * const argv[] + , int & experiment_number + , std::string & filepoints + , std::string & landmark_file + , std::string & experiment_name + , int & nbL + , double & alpha2_s + , double & alpha2_w + , double & mu_epsilon + , int & dim_max + , std::vector & desired_homology + , double & min_persistence) { + namespace po = boost::program_options; + po::options_description hidden("Hidden options"); + hidden.add_options() + ("option", po::value(& experiment_number), + "Experiment id.") + ("input-file", po::value(&filepoints), + "Name of file containing a point set. Format is one point per line: X1 ... Xd "); + + po::options_description visible("Allowed options", 100); + visible.add_options() + ("help,h", "produce help message") + ("output-file,o", po::value(&experiment_name)->default_value("witness"), + "The prefix of all the output files. Default is 'witness'") + ("landmarks,L", po::value(&nbL)->default_value(0), + "Number of landmarks.") + ( "landmark-file,l", po::value(&landmark_file), + "Name of a fike containing landmarks") + ("alpha2_s,A", po::value(&alpha2_s)->default_value(0), + "Relaxation parameter for the strong complex.") + ("alpha2_w,a", po::value(&alpha2_w)->default_value(0), + "Relaxation parameter for the weak complex.") + ("mu_epsilon,e", po::value(&mu_epsilon)->default_value(0), + "Sparsification parameter.") + ("cpx-dimension,d", po::value(&dim_max)->default_value(1), + "Maximal dimension of the Witness complex we want to compute.") + + ("homology,H", po::value>(&desired_homology)->multitoken(), + "The desired Betti numbers.") + ("min-persistence,m", po::value(&min_persistence), + "Minimal lifetime of homology feature to be recorded. Default is 0. Enter a negative value to see zero length intervals"); + + po::positional_options_description pos; + pos.add("option", 1); + pos.add("input-file", 2); + + po::options_description all; + all.add(visible).add(hidden); + + po::variables_map vm; + po::store(po::command_line_parser(argc, argv). + options(all).positional(pos).run(), vm); + po::notify(vm); + + if (vm.count("help") || !vm.count("input-file")) { + std::cout << std::endl; + std::cout << "Compute the persistent homology with coefficient field Z/3Z \n"; + std::cout << "of a Strong relaxed witness complex defined on a set of input points.\n \n"; + std::cout << "The output diagram contains one bar per line, written with the convention: \n"; + std::cout << " p dim b d \n"; + std::cout << "where dim is the dimension of the homological feature,\n"; + std::cout << "b and d are respectively the birth and death of the feature and \n"; + std::cout << "p is the characteristic of the field Z/pZ used for homology coefficients." << std::endl << std::endl; + + std::cout << "Usage: " << argv[0] << " [options] input-file" << std::endl << std::endl; + std::cout << visible << std::endl; + std::abort(); + } +} + + + + + /** * \brief Customized version of read_points * which takes into account a possible nbP first line @@ -175,6 +260,7 @@ void rw_experiment(Point_Vector & point_vector, int nbL, FT alpha2, int limD, FT Gudhi::witness_complex::Dim_lists simplices(simplex_tree, limD); // std::vector simplices; + std::cout << "Starting collapses...\n"; simplices.collapse(); simplices.output_simplices(); @@ -204,6 +290,21 @@ void rw_experiment(Point_Vector & point_vector, int nbL, FT alpha2, int limD, FT chi += 1-2*(simplex_tree.dimension(sh)%2); std::cout << "Euler characteristic is " << chi << std::endl; write_witness_mesh(point_vector, landmarks_ind, collapsed_tree, collapsed_tree.complex_simplex_range(), false, true, mesh_filename+"_after_collapse.mesh"); + Gudhi::Good_links gl(collapsed_tree); + if (gl.complex_is_pseudomanifold()) + std::cout << "Collapsed complex is a pseudomanifold.\n"; + else + std::cout << "Collapsed complex is NOT a pseudomanifold.\n"; + bool good = true; + for (auto v: collapsed_tree.complex_vertex_range()) + if (!gl.has_good_link(v)) { + std::cout << "Bad link around " << v << std::endl; + good = false; + } + if (good) + std::cout << "All links are good.\n"; + else + std::cout << "There are bad links.\n"; } void rips_experiment(Point_Vector & points, double threshold, int dim_max) @@ -267,83 +368,6 @@ int experiment0 (int argc, char * const argv[]) return 0; } -int experiment1 (int argc, char * const argv[]) -{ - if (argc != 8) { - std::cerr << "Usage: " << argv[0] - << " 1 file_name nbL alpha mu_epsilon limD experiment_name\n"; - return 0; - } - /* - boost::filesystem::path p; - for (; argc > 2; --argc, ++argv) - p /= argv[1]; - */ - - std::string file_name = argv[2]; - int nbL = atoi(argv[3]), limD = atoi(argv[6]); - double alpha2 = atof(argv[4]), mu_epsilon = atof(argv[5]); - std::string experiment_name = argv[7]; - - // 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"; - - STree simplex_tree; - std::vector > knn; - std::vector > distances; - std::vector landmarks; - Gudhi::witness_complex::landmark_choice_by_sparsification(point_vector, nbL, mu_epsilon, landmarks); - Gudhi::witness_complex::build_distance_matrix(point_vector, // aka witnesses - landmarks, // aka landmarks - alpha2, - limD, - knn, - distances); - WRWit rw(distances, - knn, - simplex_tree, - nbL, - 0, - limD); - std::vector landmarks_ind(nbL); - for (unsigned i = 0; i != distances.size(); ++i) { - if (distances[i][0] == 0) - landmarks_ind[knn[i][0]] = i; - } - - write_witness_mesh(point_vector, landmarks_ind, simplex_tree, simplex_tree.complex_simplex_range(), false, true, experiment_name+"_witness_complex.mesh"); - - int chi = 0; - for (auto sh: simplex_tree.complex_simplex_range()) - chi += 1-2*(simplex_tree.dimension(sh)%2); - std::cout << "Euler characteristic is " << chi << std::endl; - - rw_experiment(point_vector, nbL, alpha2, limD, mu_epsilon, experiment_name); - - // 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; - // } - // } - return 0; -} /******************************************************************************************** * Length of the good interval experiment @@ -380,17 +404,25 @@ double good_interval_length(const std::vector & desired_homology, STree & s std::string line; std::vector pers_endpoints; while (getline(in_stream, line)) { - int p, dim; + unsigned p, dim; double alpha_start, alpha_end; std::istringstream iss(line); iss >> p >> dim >> alpha_start >> alpha_end; + if (iss.fail()) + alpha_end = alpha2; + //std::cout << p << " " << dim << " " << alpha_start << " " << alpha_end << "\n"; + //if (dim < desired_homology.size()+1) if (alpha_start != alpha_end) { - if (alpha_end < alpha_start) - alpha_end = alpha2; + // if (alpha_end < alpha_start) + // alpha_end = alpha2; pers_endpoints.push_back(Pers_endpoint(alpha_start, true, dim)); pers_endpoints.push_back(Pers_endpoint(alpha_end, false, dim)); + std::cout << p << " " << dim << " " << alpha_start << " " << alpha_end << "\n"; } } + std::cout << "desired_homology.size() = " << desired_homology.size() << "\n"; + for (auto nd: desired_homology) + std::cout << nd << std::endl; std::cout << "Pers_endpoints.size = " << pers_endpoints.size() << std::endl; in_stream.close(); std::sort(pers_endpoints.begin(), @@ -404,13 +436,20 @@ double good_interval_length(const std::vector & desired_homology, STree & s std::cout << p.alpha << " " << p.dim << " " << p.start << "\n"; } */ - std::vector current_homology(nbL-1,0); - current_homology[0] = 1; // for the compulsary "0 0 inf" entry + std::vector current_homology(desired_homology.size(),0); + //current_homology[0] = 1; // for the compulsary "0 0 inf" entry double good_start = 0, good_end = 0; double sum_intervals = 0; int num_pieces = 0; bool interval_in_process = (desired_homology == current_homology); for (auto p: pers_endpoints) { + /* + std::cout << p.alpha << " " << p.dim << " "; + if (p.start) + std::cout << "s\n"; + else + std::cout << "e\n"; + */ /* std::cout << "Treating " << p.alpha << " " << p.dim << " " << p.start << " ["; @@ -450,66 +489,114 @@ void run_comparison(std::vector > const & knn, std::vector > const & distances, Point_Vector & points, unsigned nbL, - double alpha2, + unsigned limD, + double alpha2_s, + double alpha2_w, std::vector& desired_homology) { clock_t start, end; STree simplex_tree; - alpha2 = 0.55; - std::cout << "alpha2 = " << alpha2 << "\n"; + //std::cout << "alpha2 = " << alpha2_s << "\n"; start = clock(); SRWit srwit(distances, knn, simplex_tree, nbL, - alpha2, - nbL-1); + alpha2_s, + limD); end = clock(); std::cout << "SRWit.size = " << simplex_tree.num_simplices() << std::endl; - simplex_tree.set_dimension(nbL-1); + simplex_tree.set_dimension(desired_homology.size()); std::cout << "Good homology interval length for SRWit is " - << good_interval_length(desired_homology, simplex_tree, alpha2) << "\n"; + << good_interval_length(desired_homology, simplex_tree, alpha2_s) << "\n"; std::cout << "Time: " << static_cast(end - start) / CLOCKS_PER_SEC << " s. \n"; + int chi = 0; + for (auto sh: simplex_tree.complex_simplex_range()) + chi += 1-2*(simplex_tree.dimension(sh)%2); + std::cout << "Euler characteristic is " << chi << std::endl; + STree simplex_tree2; - alpha2 = 0.35; - std::cout << "alpha2 = " << alpha2 << "\n"; + std::cout << "alpha2 = " << alpha2_w << "\n"; start = clock(); WRWit wrwit(distances, knn, simplex_tree2, nbL, - alpha2, - nbL-1); + alpha2_w, + limD); end = clock(); std::cout << "WRWit.size = " << simplex_tree2.num_simplices() << std::endl; - simplex_tree.set_dimension(nbL-1); + simplex_tree2.set_dimension(nbL-1); std::cout << "Good homology interval length for WRWit is " - << good_interval_length(desired_homology, simplex_tree2, alpha2) << "\n"; + << good_interval_length(desired_homology, simplex_tree2, alpha2_w) << "\n"; std::cout << "Time: " << static_cast(end - start) / CLOCKS_PER_SEC << " s. \n"; + chi = 0; + for (auto sh: simplex_tree2.complex_simplex_range()) + chi += 1-2*(simplex_tree2.dimension(sh)%2); + std::cout << "Euler characteristic is " << chi << std::endl; + + //write_witness_mesh(points, landmarks_ind, simplex_tree2, simplex_tree2.complex_simplex_range(), false, true, "wrwit.mesh"); - STree simplex_tree3; - std::cout << "Enter alpha2 for Rips: "; - std::cin >> alpha2; - std::cout << "\n"; - std::cout << "points.size() = " << points.size() << "\n"; - start = clock(); - rips(points, - 0.2, - nbL-1, - simplex_tree3); - end = clock(); - std::cout << "Rips.size = " << simplex_tree3.num_simplices() << std::endl; - simplex_tree.set_dimension(nbL-1); - std::cout << "Good homology interval length for WRWit is " - << good_interval_length(desired_homology, simplex_tree3, alpha2) << "\n"; - std::cout << "Time: " << static_cast(end - start) / CLOCKS_PER_SEC << " s. \n"; } +int experiment1 (int argc, char * const argv[]) +{ + /* + boost::filesystem::path p; + for (; argc > 2; --argc, ++argv) + p /= argv[1]; + */ + + // std::string file_name = argv[2]; + // int nbL = atoi(argv[3]), limD = atoi(argv[6]); + // double alpha2 = atof(argv[4]), mu_epsilon = atof(argv[5]); + // std::string experiment_name = argv[7]; + + int option = 1; + std::string file_name, landmark_file; + int nbL = 0, limD; + double alpha2_s, alpha2_w, mu_epsilon, min_pers; + std::string experiment_name; + std::vector desired_homology = {1}; + std::vector landmarks; + + program_options(argc, argv, option, file_name, landmark_file, experiment_name, nbL, alpha2_s, alpha2_w, mu_epsilon, limD, desired_homology, min_pers); + + // 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"; + //std::cout << "Limit dimension for the complexes is " << limD << ".\n"; + + if (landmark_file == "") + Gudhi::witness_complex::landmark_choice_by_sparsification(point_vector, nbL, mu_epsilon, landmarks); + //Gudhi::witness_complex::landmark_choice_by_random_knn(point_vector, nbL, alpha2_s, limD, knn, distances); + else + read_points_cust(landmark_file, landmarks); + nbL = landmarks.size(); + STree simplex_tree; + std::vector > knn; + std::vector > distances; + + //Gudhi::witness_complex::landmark_choice_by_sparsification(point_vector, nbL, mu_epsilon, landmarks); + Gudhi::witness_complex::build_distance_matrix(point_vector, // aka witnesses + landmarks, // aka landmarks + alpha2_s, + limD, + knn, + distances); + + run_comparison(knn, distances, point_vector, nbL, limD, alpha2_s, alpha2_w, desired_homology); + return 0; +} + + int experiment2(int argc, char * const argv[]) { for (unsigned d = 3; d < 4; d++) { @@ -524,7 +611,6 @@ int experiment2(int argc, char * const argv[]) case 4: alpha2 = 1.4; break; default: alpha2 = 1.4; break; } - std::cout << "alpha2 = " << alpha2 << "\n"; unsigned nbL = 20; std::vector desired_homology(nbL-1,0); desired_homology[0] = 1; desired_homology[d] = 1; @@ -536,19 +622,22 @@ int experiment2(int argc, char * const argv[]) std::cout << "Running test S"<< d <<", |W|=" << nbW << ", |L|=" << nbL << std::endl; generate_points_sphere(point_vector, i*1000, d+1); std::vector landmarks; - Gudhi::witness_complex::landmark_choice_by_sparsification(point_vector, nbL, mu_epsilon, landmarks); + Gudhi::witness_complex::landmark_choice_by_sparsification(point_vector, nbL, mu_epsilon, landmarks); + std::vector > knn; std::vector > distances; - + + std::cout << "|L| after sparsification: " << landmarks.size() << "\n"; + Gudhi::witness_complex::build_distance_matrix(point_vector, // aka witnesses landmarks, // aka landmarks alpha2, nbL-1, knn, distances); - run_comparison(knn, distances, point_vector, nbL, alpha2, desired_homology); + run_comparison(knn, distances, point_vector, nbL, nbL-1, alpha2, alpha2, desired_homology); } } /* @@ -588,7 +677,8 @@ int experiment2(int argc, char * const argv[]) int experiment3(int argc, char * const argv[]) { - // COLLAPSES EXPERIMENT + // Both witnesses and landmarks are given as input + return 0; } @@ -609,6 +699,9 @@ int main (int argc, char * const argv[]) case 2 : return experiment2(argc, argv); break; + case 3 : + return experiment3(argc, argv); + break; default : output_experiment_information(argv[0]); return 1; diff --git a/src/Witness_complex/example/generators.h b/src/Witness_complex/example/generators.h index 23915546..731a52b0 100644 --- a/src/Witness_complex/example/generators.h +++ b/src/Witness_complex/example/generators.h @@ -25,10 +25,12 @@ #include #include +#include #include #include #include +#include typedef CGAL::Epick_d K; typedef K::FT FT; @@ -144,12 +146,21 @@ void generate_points_sphere(Point_Vector& W, int nbP, int dim) { W.push_back(*rp++); } -/** \brief Generate nbP points on a d-torus +/** \brief Generate nbP points on a (flat) d-torus embedded in R^{2d} * */ void generate_points_torus(Point_Vector& W, int nbP, int dim) { - CGAL::Random_points_on_sphere_d rp(2, 1); - + CGAL::Random rand; + const double pi = std::acos(-1); + for (int i = 0; i < nbP; i++) { + std::vector point; + for (int j = 0; j < dim; j++) { + double alpha = rand.uniform_real((double)0, 2*pi); + point.push_back(sin(alpha)); + point.push_back(cos(alpha)); + } + W.push_back(Point_d(point)); + } } #endif // EXAMPLE_WITNESS_COMPLEX_GENERATORS_H_ diff --git a/src/Witness_complex/example/output.h b/src/Witness_complex/example/output.h index 0e74387a..d3f534af 100644 --- a/src/Witness_complex/example/output.h +++ b/src/Witness_complex/example/output.h @@ -73,9 +73,12 @@ void write_edges(std::string file_name, STree& witness_complex, Point_Vector& la ofs.close(); } -/** \brief Write triangles (tetrahedra in 3d) of a witness - * complex in a file, compatible with medit. - * l_is_v = landmark is vertex +/** \brief Write triangles (tetrahedra in 3d) of a simplicial complex in a file, compatible with medit. + * `landmarks_ind` represents the set of landmark indices in W + * `st` is the Simplex_tree to be visualized, + * `shr` is the Simplex_handle_range of simplices in `st` to be visualized + * `is2d` should be true if the simplicial complex is 2d, false if 3d + * `l_is_v` = landmark is vertex */ template diff --git a/src/Witness_complex/include/gudhi/Good_links.h b/src/Witness_complex/include/gudhi/Good_links.h index a011c032..c5e993e5 100644 --- a/src/Witness_complex/include/gudhi/Good_links.h +++ b/src/Witness_complex/include/gudhi/Good_links.h @@ -50,23 +50,400 @@ #include namespace Gudhi { + +template < typename Simplicial_complex > +class Good_links { + + typedef typename Simplicial_complex::Simplex_handle Simplex_handle; + typedef typename Simplicial_complex::Vertex_handle Vertex_handle; + typedef std::vector Vertex_vector; + + // graph is an adjacency list + typedef typename boost::adjacency_list Adj_graph; + // map that gives to a certain simplex its node in graph and its dimension + //typedef std::pair Reference; + typedef boost::graph_traits::vertex_descriptor Vertex_t; + typedef boost::graph_traits::edge_descriptor Edge_t; + typedef boost::graph_traits::adjacency_iterator Adj_it; + typedef std::pair Out_edge_it; + + typedef boost::container::flat_map Graph_map; + typedef boost::container::flat_map Inv_graph_map; + +public: + Good_links(Simplicial_complex& sc): sc_(sc) + { + int dim = 0; + for (auto sh: sc_.complex_simplex_range()) + if (sc_.dimension(sh) > dim) + dim = sc_.dimension(sh); + dimension = dim; + count_good = Vertex_vector(dim); + count_bad = Vertex_vector(dim); + } + +private: + + Simplicial_complex& sc_; + unsigned dimension; + std::vector count_good; + std::vector count_bad; + + void add_vertices_to_link_graph(Vertex_vector& star_vertices, + typename Vertex_vector::iterator curr_v, + Adj_graph& adj_graph, + Graph_map& d_map, + Graph_map& f_map, + Vertex_vector& curr_simplex, + int curr_d, + int link_dimension) + { + Simplex_handle sh; + Vertex_t vert; + typename Vertex_vector::iterator it; + //std::pair resPair; + //typename Graph_map::iterator resPair; + //Add vertices + //std::cout << "Entered add vertices\n"; + for (it = curr_v; it != star_vertices.end(); ++it) { + curr_simplex.push_back(*it); //push next vertex in question + curr_simplex.push_back(star_vertices[0]); //push the center of the star + /* + std::cout << "Searching for "; + print_vector(curr_simplex); + std::cout << " curr_dim " << curr_d << " d " << dimension << ""; + */ + Vertex_vector curr_simplex_copy(curr_simplex); + sh = sc_.find(curr_simplex_copy); //a simplex of the star + curr_simplex.pop_back(); //pop the center of the star + curr_simplex_copy = Vertex_vector(curr_simplex); + if (sh != sc_.null_simplex()) { + //std::cout << " added\n"; + if (curr_d == link_dimension) { + sh = sc_.find(curr_simplex_copy); //a simplex of the link + assert(sh != sc_.null_simplex()); //ASSERT! + vert = boost::add_vertex(adj_graph); + d_map.emplace(sh,vert); + } + else { + if (curr_d == link_dimension-1) { + sh = sc_.find(curr_simplex_copy); //a simplex of the link + assert(sh != sc_.null_simplex()); + vert = boost::add_vertex(adj_graph); + f_map.emplace(sh,vert); + } + + //delete (&curr_simplex_copy); //Just so it doesn't stack + add_vertices_to_link_graph(star_vertices, + it+1, + adj_graph, + d_map, + f_map, + curr_simplex, + curr_d+1, link_dimension); + } + } + /* + else + std::cout << "\n"; + */ + curr_simplex.pop_back(); //pop the vertex in question + } + } -namespace witness_complex { + void add_edges_to_link_graph(Adj_graph& adj_graph, + Graph_map& d_map, + Graph_map& f_map) + { + Simplex_handle sh; + // Add edges + //std::cout << "Entered add edges:\n"; + typename Graph_map::iterator map_it; + for (auto d_map_pair : d_map) { + //std::cout << "*"; + sh = d_map_pair.first; + Vertex_t d_vert = d_map_pair.second; + for (auto facet_sh : sc_.boundary_simplex_range(sh)) + //for (auto f_map_it : f_map) + { + //std::cout << "'"; + map_it = f_map.find(facet_sh); + //We must have all the facets in the graph at this point + assert(map_it != f_map.end()); + Vertex_t f_vert = map_it->second; + //std::cout << "Added edge " << sh->first << "-" << map_it->first->first << "\n"; + boost::add_edge(d_vert,f_vert,adj_graph); + } + } + } - /** \addtogroup simplex_tree - * Witness complex is a simplicial complex defined on two sets of points in \f$\mathbf{R}^D\f$: - * \f$W\f$ set of witnesses and \f$L \subseteq W\f$ set of landmarks. The simplices are based on points in \f$L\f$ - * and a simplex belongs to the witness complex if and only if it is witnessed (there exists a point \f$w \in W\f$ such that - * w is closer to the vertices of this simplex than others) and all of its faces are witnessed as well. + + + /* \brief Verifies if the simplices formed by vertices given by link_vertices + * form a pseudomanifold. + * The idea is to make a bipartite graph, where vertices are the d- and (d-1)-dimensional + * faces and edges represent adjacency between them. + */ + bool link_is_pseudomanifold(Vertex_vector& star_vertices, + int dimension) + { + Adj_graph adj_graph; + Graph_map d_map, f_map; + // d_map = map for d-dimensional simplices + // f_map = map for its facets + Vertex_vector init_vector = {}; + add_vertices_to_link_graph(star_vertices, + star_vertices.begin()+1, + adj_graph, + d_map, + f_map, + init_vector, + 0, dimension); + //std::cout << "DMAP_SIZE: " << d_map.size() << "\n"; + //std::cout << "FMAP_SIZE: " << f_map.size() << "\n"; + add_edges_to_link_graph(adj_graph, + d_map, + f_map); + for (auto f_map_it : f_map) { + //std::cout << "Degree of " << f_map_it.first->first << " is " << boost::out_degree(f_map_it.second, adj_graph) << "\n"; + if (boost::out_degree(f_map_it.second, adj_graph) != 2){ + count_bad[dimension]++; + return false; + } + } + // At this point I know that all (d-1)-simplices are adjacent to exactly 2 d-simplices + // What is left is to check the connexity + //std::vector components(boost::num_vertices(adj_graph)); + return true; //Forget the connexity + //return (boost::connected_components(adj_graph, &components[0]) == 1); + } + + int star_dim(Vertex_vector& star_vertices, + typename Vertex_vector::iterator curr_v, + int curr_d, + Vertex_vector& curr_simplex, + typename std::vector< int >::iterator curr_dc) + { + //std::cout << "Entered star_dim for " << *(curr_v-1) << "\n"; + Simplex_handle sh; + int final_d = curr_d; + typename Vertex_vector::iterator it; + typename Vertex_vector::iterator dc_it; + //std::cout << "Current vertex is " << + for (it = curr_v, dc_it = curr_dc; it != star_vertices.end(); ++it, ++dc_it) + { + curr_simplex.push_back(*it); + Vertex_vector curr_simplex_copy(curr_simplex); + /* + std::cout << "Searching for "; + print_vector(curr_simplex); + std::cout << " curr_dim " << curr_d << " final_dim " << final_d; + */ + sh = sc_.find(curr_simplex_copy); //Need a copy because find sorts the vector and I want star center to be the first + if (sh != sc_.null_simplex()) + { + //std::cout << " -> " << *it << "\n"; + int d = star_dim(star_vertices, + it+1, + curr_d+1, + curr_simplex, + dc_it); + if (d >= final_d) + { + final_d = d; + //std::cout << d << " "; + //print_vector(curr_simplex); + //std::cout << std::endl; + } + if (d >= *dc_it) + *dc_it = d; + } + /* + else + std::cout << "\n"; + */ + curr_simplex.pop_back(); + } + return final_d; + } + +public: + + /** \brief Returns true if the link is good */ -template< class Simplicial_complex > -bool good_links(Simplicial_complex& sc_) -{ + bool has_good_link(Vertex_handle v) + { + typedef Vertex_vector typeVectorVertex; + typeVectorVertex star_vertices; + // Fill star_vertices + star_vertices.push_back(v); + for (auto u: sc_.complex_vertex_range()) + { + typeVectorVertex edge = {u,v}; + if (u != v && sc_.find(edge) != sc_.null_simplex()) + star_vertices.push_back(u); + } + // Find the dimension + typeVectorVertex init_simplex = {star_vertices[0]}; + bool is_pure = true; + std::vector dim_coface(star_vertices.size(), 1); + int d = star_dim(star_vertices, + star_vertices.begin()+1, + 0, + init_simplex, + dim_coface.begin()+1) - 1; //link_dim = star_dim - 1 + + assert(init_simplex.size() == 1); + // if (!is_pure) + // std::cout << "Found an impure star around " << v << "\n"; + for (int dc: dim_coface) + is_pure = (dc == dim_coface[0]); + /* + if (d == count_good.size()) + { + std::cout << "Found a star of dimension " << (d+1) << " around " << v << "\nThe star is "; + print_vector(star_vertices); std::cout << std::endl; + } + */ + //if (d == -1) count_bad[0]++; + bool b= (is_pure && link_is_pseudomanifold(star_vertices,d)); + if (d != -1) {if (b) count_good[d]++; else count_bad[d]++;} + if (!is_pure) count_bad[0]++; + return (d != -1 && b && is_pure); + } + +private: + void add_max_simplices_to_graph(Vertex_vector& star_vertices, + typename Vertex_vector::iterator curr_v, + Adj_graph& adj_graph, + Graph_map& d_map, + Graph_map& f_map, + Inv_graph_map& inv_d_map, + Vertex_vector& curr_simplex, + int curr_d, + int link_dimension) + { + Simplex_handle sh; + Vertex_t vert; + typename Vertex_vector::iterator it; + //std::pair resPair; + //typename Graph_map::iterator resPair; + //Add vertices + //std::cout << "Entered add vertices\n"; + for (it = curr_v; it != star_vertices.end(); ++it) { + curr_simplex.push_back(*it); //push next vertex in question + //curr_simplex.push_back(star_vertices[0]); //push the center of the star + /* + std::cout << "Searching for "; + print_vector(curr_simplex); + std::cout << " curr_dim " << curr_d << " d " << dimension << ""; + */ + Vertex_vector curr_simplex_copy(curr_simplex); + sh = sc_.find(curr_simplex_copy); //a simplex of the star + //curr_simplex.pop_back(); //pop the center of the star + curr_simplex_copy = Vertex_vector(curr_simplex); + if (sh != sc_.null_simplex()) { + //std::cout << " added\n"; + if (curr_d == link_dimension) { + sh = sc_.find(curr_simplex_copy); //a simplex of the link + assert(sh != sc_.null_simplex()); //ASSERT! + vert = boost::add_vertex(adj_graph); + d_map.emplace(sh,vert); + inv_d_map.emplace(vert,sh); + } + else { + if (curr_d == link_dimension-1) { + sh = sc_.find(curr_simplex_copy); //a simplex of the link + assert(sh != sc_.null_simplex()); + vert = boost::add_vertex(adj_graph); + f_map.emplace(sh,vert); + } + //delete (&curr_simplex_copy); //Just so it doesn't stack + add_max_simplices_to_graph(star_vertices, + it+1, + adj_graph, + d_map, + f_map, + inv_d_map, + curr_simplex, + curr_d+1, + link_dimension); + } + } + /* + else + std::cout << "\n"; + */ + curr_simplex.pop_back(); //pop the vertex in question + } + } + +public: + bool complex_is_pseudomanifold() + { + Adj_graph adj_graph; + Graph_map d_map, f_map; + // d_map = map for d-dimensional simplices + // f_map = map for its facets + Inv_graph_map inv_d_map; + Vertex_vector init_vector = {}; + std::vector star_vertices; + for (int v: sc_.complex_vertex_range()) + star_vertices.push_back(v); + add_max_simplices_to_graph(star_vertices, + star_vertices.begin(), + adj_graph, + d_map, + f_map, + inv_d_map, + init_vector, + 0, dimension); + //std::cout << "DMAP_SIZE: " << d_map.size() << "\n"; + //std::cout << "FMAP_SIZE: " << f_map.size() << "\n"; + add_edges_to_link_graph(adj_graph, + d_map, + f_map); + for (auto f_map_it : f_map) { + //std::cout << "Degree of " << f_map_it.first->first << " is " << boost::out_degree(f_map_it.second, adj_graph) << "\n"; + if (boost::out_degree(f_map_it.second, adj_graph) != 2) { + // if (boost::out_degree(f_map_it.second, adj_graph) >= 3) { + // // std::cout << "This simplex has 3+ cofaces: "; + // // for(auto v : sc_.simplex_vertex_range(f_map_it.first)) + // // std::cout << v << " "; + // // std::cout << std::endl; + // Adj_it ai, ai_end; + // for (std::tie(ai, ai_end) = boost::adjacent_vertices(f_map_it.second, adj_graph); ai != ai_end; ++ai) { + // auto it = inv_d_map.find(*ai); + // assert (it != inv_d_map.end()); + // typename Simplicial_complex::Simplex_handle sh = it->second; + // for(auto v : sc_.simplex_vertex_range(sh)) + // std::cout << v << " "; + // std::cout << std::endl; + // } + // } + count_bad[dimension]++; + return false; + } + } + // At this point I know that all (d-1)-simplices are adjacent to exactly 2 d-simplices + // What is left is to check the connexity + //std::vector components(boost::num_vertices(adj_graph)); + return true; //Forget the connexity + //return (boost::connected_components(adj_graph, &components[0]) == 1); + } -} + int number_good_links(int dim) + { + return count_good[dim]; + } -} // namespace witness_complex + int number_bad_links(int dim) + { + return count_bad[dim]; + } +}; + } // namespace Gudhi #endif -- cgit v1.2.3