From 1bcd448f6e217655fae4bfbf62d8d3b6caa52503 Mon Sep 17 00:00:00 2001 From: skachano Date: Fri, 7 Oct 2016 09:49:47 +0000 Subject: Added the possibility to limit dimension git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/relaxed-witness@1676 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: 16ec8706ffb981ff892960468cc4b809e61612e2 --- .../example/example_strong_witness_persistence.cpp | 10 +++++++--- .../example/example_witness_complex_off.cpp | 6 +++--- .../example_witness_complex_persistence.cpp | 12 +++++++---- .../include/gudhi/Active_witness/Active_witness.h | 3 +++ .../gudhi/Active_witness/Active_witness_iterator.h | 7 +++++++ .../include/gudhi/Strong_witness_complex.h | 23 ++++++++++++++-------- .../include/gudhi/Witness_complex.h | 23 +++++++++++++--------- 7 files changed, 57 insertions(+), 27 deletions(-) (limited to 'src/Witness_complex') diff --git a/src/Witness_complex/example/example_strong_witness_persistence.cpp b/src/Witness_complex/example/example_strong_witness_persistence.cpp index bf6b4011..4e06867d 100644 --- a/src/Witness_complex/example/example_strong_witness_persistence.cpp +++ b/src/Witness_complex/example/example_strong_witness_persistence.cpp @@ -53,17 +53,18 @@ void program_options(int argc, char * argv[] , std::string & filediag , Filtration_value & max_squared_alpha , int & p + , int & dim_max , Filtration_value & min_persistence); int main(int argc, char * argv[]) { std::string file_name; std::string filediag; Filtration_value max_squared_alpha; - int p, nbL; + int p, nbL, lim_d; Filtration_value min_persistence; SimplexTree simplex_tree; - program_options(argc, argv, nbL, file_name, filediag, max_squared_alpha, p, min_persistence); + program_options(argc, argv, nbL, file_name, filediag, max_squared_alpha, p, lim_d, min_persistence); // Extract the points from the file file_name Point_vector witnesses, landmarks; @@ -119,6 +120,7 @@ void program_options(int argc, char * argv[] , std::string & filediag , Filtration_value & max_squared_alpha , int & p + , int & dim_max , Filtration_value & min_persistence) { namespace po = boost::program_options; @@ -141,7 +143,9 @@ void program_options(int argc, char * argv[] ("field-charac,p", po::value(&p)->default_value(11), "Characteristic p of the coefficient field Z/pZ for computing homology.") ("min-persistence,m", po::value(&min_persistence)->default_value(0), - "Minimal lifetime of homology feature to be recorded. Default is 0. Enter a negative value to see zero length intervals"); + "Minimal lifetime of homology feature to be recorded. Default is 0. Enter a negative value to see zero length intervals") + ("cpx-dimension,d", po::value(&dim_max)->default_value(std::numeric_limits::max()), + "Maximal dimension of the strong witness complex we want to compute."); po::positional_options_description pos; pos.add("input-file", 1); diff --git a/src/Witness_complex/example/example_witness_complex_off.cpp b/src/Witness_complex/example/example_witness_complex_off.cpp index 849a0a40..37646faf 100644 --- a/src/Witness_complex/example/example_witness_complex_off.cpp +++ b/src/Witness_complex/example/example_witness_complex_off.cpp @@ -45,12 +45,12 @@ typedef std::vector< Point_d > Point_vector; int main(int argc, char * const argv[]) { if (argc != 4) { std::cerr << "Usage: " << argv[0] - << " path_to_point_file nbL alpha^2\n"; + << " path_to_point_file number_of_landmarks max_squared_alpha limit_dimension\n"; return 0; } std::string file_name = argv[1]; - int nbL = atoi(argv[2]); + int nbL = atoi(argv[2]), lim_dim = atoi(argv[4]); double alpha2 = atof(argv[3]); clock_t start, end; Gudhi::Simplex_tree<> simplex_tree; @@ -77,7 +77,7 @@ int main(int argc, char * const argv[]) { point_vector.begin(), point_vector.end()); - witness_complex.create_complex(simplex_tree, alpha2); + witness_complex.create_complex(simplex_tree, alpha2, lim_dim); end = clock(); std::cout << "Witness complex took " << static_cast(end - start) / CLOCKS_PER_SEC << " s. \n"; diff --git a/src/Witness_complex/example/example_witness_complex_persistence.cpp b/src/Witness_complex/example/example_witness_complex_persistence.cpp index 9b06b504..a4388047 100644 --- a/src/Witness_complex/example/example_witness_complex_persistence.cpp +++ b/src/Witness_complex/example/example_witness_complex_persistence.cpp @@ -53,17 +53,18 @@ void program_options(int argc, char * argv[] , std::string & filediag , Filtration_value & max_squared_alpha , int & p + , int & dim_max , Filtration_value & min_persistence); int main(int argc, char * argv[]) { std::string file_name; std::string filediag; Filtration_value max_squared_alpha; - int p, nbL; + int p, nbL, lim_d; Filtration_value min_persistence; SimplexTree simplex_tree; - program_options(argc, argv, nbL, file_name, filediag, max_squared_alpha, p, min_persistence); + program_options(argc, argv, nbL, file_name, filediag, max_squared_alpha, p, lim_d, min_persistence); // Extract the points from the file file_name Point_vector witnesses, landmarks; @@ -119,6 +120,7 @@ void program_options(int argc, char * argv[] , std::string & filediag , Filtration_value & max_squared_alpha , int & p + , int & dim_max , Filtration_value & min_persistence) { namespace po = boost::program_options; @@ -141,8 +143,10 @@ void program_options(int argc, char * argv[] ("field-charac,p", po::value(&p)->default_value(11), "Characteristic p of the coefficient field Z/pZ for computing homology.") ("min-persistence,m", po::value(&min_persistence)->default_value(0), - "Minimal lifetime of homology feature to be recorded. Default is 0. Enter a negative value to see zero length intervals"); - + "Minimal lifetime of homology feature to be recorded. Default is 0. Enter a negative value to see zero length intervals") + ("cpx-dimension,d", po::value(&dim_max)->default_value(std::numeric_limits::max()), + "Maximal dimension of the weak witness complex we want to compute."); + po::positional_options_description pos; pos.add("input-file", 1); diff --git a/src/Witness_complex/include/gudhi/Active_witness/Active_witness.h b/src/Witness_complex/include/gudhi/Active_witness/Active_witness.h index 9ae41a69..e52410e4 100644 --- a/src/Witness_complex/include/gudhi/Active_witness/Active_witness.h +++ b/src/Witness_complex/include/gudhi/Active_witness/Active_witness.h @@ -31,6 +31,9 @@ namespace Gudhi { namespace witness_complex { + /** \brief Class representing a list of nearest neighbors to a given witness. + * \detail Every element is a pair of a landmark identifier and the squared distance to it. + */ template< typename Id_distance_pair, typename INS_range > class Active_witness { diff --git a/src/Witness_complex/include/gudhi/Active_witness/Active_witness_iterator.h b/src/Witness_complex/include/gudhi/Active_witness/Active_witness_iterator.h index 5d4f3d75..9c96f7e8 100644 --- a/src/Witness_complex/include/gudhi/Active_witness/Active_witness_iterator.h +++ b/src/Witness_complex/include/gudhi/Active_witness/Active_witness_iterator.h @@ -31,6 +31,13 @@ namespace Gudhi { namespace witness_complex { + /** \brief Iterator in the nearest landmark list. + * \detail After the iterator reaches the end of the list, + * the list is augmented by a (nearest landmark, distance) pair if possible. + * If all the landmarks are present in the list, iterator returns the specific end value + * of the corresponding 'Active_witness' object. + */ + template< typename Active_witness, typename Id_distance_pair, typename INS_iterator > diff --git a/src/Witness_complex/include/gudhi/Strong_witness_complex.h b/src/Witness_complex/include/gudhi/Strong_witness_complex.h index 3a862d71..e125d307 100644 --- a/src/Witness_complex/include/gudhi/Strong_witness_complex.h +++ b/src/Witness_complex/include/gudhi/Strong_witness_complex.h @@ -138,19 +138,26 @@ private: /** \brief Outputs the (weak) witness complex with * squared relaxation parameter 'max_alpha_square' * to simplicial complex 'complex'. + * The parameter 'limit_dimension' represents the maximal dimension of the simplicial complex + * (default value = no limit). */ template < typename SimplicialComplexForWitness > bool create_complex(SimplicialComplexForWitness& complex, - FT max_alpha_square) + FT max_alpha_square, + Landmark_id limit_dimension = std::numeric_limits::max()) { - unsigned nbL = landmarks_.size(); - unsigned complex_dim = 0; + std::size_t nbL = landmarks_.size(); + Landmark_id complex_dim = 0; if (complex.num_vertices() > 0) { - std::cerr << "Witness complex cannot create complex - complex is not empty.\n"; + std::cerr << "Strong witness complex cannot create complex - complex is not empty.\n"; return false; } if (max_alpha_square < 0) { - std::cerr << "Witness complex cannot create complex - squared relaxation parameter must be non-negative.\n"; + std::cerr << "Strong witness complex cannot create complex - squared relaxation parameter must be non-negative.\n"; + return false; + } + if (limit_dimension < 0) { + std::cerr << "Strong witness complex cannot create complex - limit dimension must be non-negative.\n"; return false; } typeVectorVertex vv; @@ -165,13 +172,13 @@ private: ActiveWitness aw(landmark_tree_.query_incremental_nearest_neighbors(w)); typeVectorVertex simplex; typename ActiveWitness::iterator aw_it = aw.begin(); - float lim_d2 = aw.begin()->second + max_alpha_square; - while (aw_it != aw.end() && aw_it->second < lim_d2) { + float lim_dist2 = aw.begin()->second + max_alpha_square; + while ((Landmark_id)simplex.size() <= limit_dimension + 1 && aw_it != aw.end() && aw_it->second < lim_dist2) { simplex.push_back(aw_it->first); complex.insert_simplex_and_subfaces(simplex, aw_it->second - aw.begin()->second); aw_it++; } - if (simplex.size() - 1 > complex_dim) + if ((Landmark_id)simplex.size() - 1 > complex_dim) complex_dim = simplex.size() - 1; } complex.set_dimension(complex_dim); diff --git a/src/Witness_complex/include/gudhi/Witness_complex.h b/src/Witness_complex/include/gudhi/Witness_complex.h index 7d8fce86..2a89306d 100644 --- a/src/Witness_complex/include/gudhi/Witness_complex.h +++ b/src/Witness_complex/include/gudhi/Witness_complex.h @@ -129,7 +129,8 @@ private: */ template < typename SimplicialComplexForWitness > bool create_complex(SimplicialComplexForWitness& complex, - FT max_alpha_square) + FT max_alpha_square, + Landmark_id limit_dimension = std::numeric_limits::max()) { std::size_t nbL = landmarks_.size(); if (complex.num_vertices() > 0) { @@ -140,6 +141,10 @@ private: std::cerr << "Witness complex cannot create complex - squared relaxation parameter must be non-negative.\n"; return false; } + if (limit_dimension < 0) { + std::cerr << "Witness complex cannot create complex - limit dimension must be non-negative.\n"; + return false; + } typeVectorVertex vv; ActiveWitnessList active_witnesses;// = new ActiveWitnessList(); for (unsigned i = 0; i != nbL; ++i) { @@ -150,13 +155,13 @@ private: complex.insert_simplex(vv, Filtration_value(0.0)); /* TODO Error if not inserted : normally no need here though*/ } - unsigned k = 1; /* current dimension in iterative construction */ + Landmark_id k = 1; /* current dimension in iterative construction */ for (auto w: witnesses_) active_witnesses.push_back(ActiveWitness(landmark_tree_.query_incremental_nearest_neighbors(w))); ActiveWitness aw_copy(active_witnesses.front()); - while (!active_witnesses.empty() && k < nbL ) { + while (!active_witnesses.empty() && k <= limit_dimension ) { typename ActiveWitnessList::iterator aw_it = active_witnesses.begin(); - std::vector simplex; + std::vector simplex; simplex.reserve(k+1); while (aw_it != active_witnesses.end()) { bool ok = add_all_faces_of_dimension(k, @@ -191,7 +196,7 @@ private: double alpha2, double norelax_dist2, typename ActiveWitness::iterator curr_l, - std::vector& simplex, + std::vector& simplex, SimplicialComplexForWitness& sc, typename ActiveWitness::iterator end) { @@ -251,17 +256,17 @@ private: * inserted_vertex is the handle of the (k+1)-th vertex witnessed by witness_id */ template < typename SimplicialComplexForWitness > - bool all_faces_in(std::vector& simplex, + bool all_faces_in(typeVectorVertex& simplex, double* filtration_value, SimplicialComplexForWitness& sc) { typedef typename SimplicialComplexForWitness::Simplex_handle Simplex_handle; - std::vector< int > facet; - for (std::vector::iterator not_it = simplex.begin(); not_it != simplex.end(); ++not_it) + typeVectorVertex facet; + for (typename typeVectorVertex::iterator not_it = simplex.begin(); not_it != simplex.end(); ++not_it) { facet.clear(); - for (std::vector::iterator it = simplex.begin(); it != simplex.end(); ++it) + for (typename typeVectorVertex::iterator it = simplex.begin(); it != simplex.end(); ++it) if (it != not_it) facet.push_back(*it); Simplex_handle facet_sh = sc.find(facet); -- cgit v1.2.3