From 5e3871c2d99249bd7560faa4e1229d3b3751b0ff Mon Sep 17 00:00:00 2001 From: ROUVREAU Vincent Date: Mon, 14 Sep 2020 11:17:40 +0200 Subject: alpha_complex_persistence utils can use weighted points. Update Alpha_kernel_d_unit_test to work --- .../utilities/alpha_complex_persistence.cpp | 101 ++++++++++++++++----- 1 file changed, 76 insertions(+), 25 deletions(-) (limited to 'src/Alpha_complex/utilities/alpha_complex_persistence.cpp') diff --git a/src/Alpha_complex/utilities/alpha_complex_persistence.cpp b/src/Alpha_complex/utilities/alpha_complex_persistence.cpp index e17831d9..e86b34e2 100644 --- a/src/Alpha_complex/utilities/alpha_complex_persistence.cpp +++ b/src/Alpha_complex/utilities/alpha_complex_persistence.cpp @@ -17,19 +17,82 @@ #include // to construct a simplex_tree from alpha complex #include +#include #include #include #include // for numeric_limits +#include +#include using Simplex_tree = Gudhi::Simplex_tree<>; using Filtration_value = Simplex_tree::Filtration_value; void program_options(int argc, char *argv[], std::string &off_file_points, bool &exact, bool &fast, - std::string &output_file_diag, Filtration_value &alpha_square_max_value, + std::string &weight_file, std::string &output_file_diag, Filtration_value &alpha_square_max_value, int &coeff_field_characteristic, Filtration_value &min_persistence); +template +std::vector read_off(const std::string &off_file_points) { + Gudhi::Points_off_reader off_reader(off_file_points); + if (!off_reader.is_valid()) { + std::cerr << "Alpha_complex - Unable to read file " << off_file_points << "\n"; + exit(-1); // ----- >> + } + return off_reader.get_point_cloud(); +} + +std::vector read_weight_file(const std::string &weight_file) { + std::vector weights; + // Read weights information from file + std::ifstream weights_ifstr(weight_file); + if (weights_ifstr.good()) { + double weight = 0.0; + // Attempt read the weight in a double format, return false if it fails + while (weights_ifstr >> weight) { + weights.push_back(weight); + } + } else { + std::cerr << "Unable to read weights file " << weight_file << std::endl; + exit(-1); + } + return weights; +} + +template +Simplex_tree create_simplex_tree(const std::string &off_file_points, const std::string &weight_file, + bool exact_version, Filtration_value alpha_square_max_value) { + Simplex_tree stree; + auto points = read_off(off_file_points); + + if (weight_file != std::string()) { + std::vector weights = read_weight_file(weight_file); + if (points.size() != weights.size()) { + std::cerr << "Alpha_complex - Inconsistency between number of points (" << points.size() + << ") and number of weights (" << weights.size() << ")" << "\n"; + exit(-1); // ----- >> + } + // Init of an alpha complex from an OFF file + Gudhi::alpha_complex::Alpha_complex alpha_complex_from_file(points, weights); + + if (!alpha_complex_from_file.create_complex(stree, alpha_square_max_value, exact_version)) { + std::cerr << "Alpha complex simplicial complex creation failed." << std::endl; + exit(-1); + } + } else { + // Init of an alpha complex from an OFF file + Gudhi::alpha_complex::Alpha_complex alpha_complex_from_file(points); + + if (!alpha_complex_from_file.create_complex(stree, alpha_square_max_value, exact_version)) { + std::cerr << "Alpha complex simplicial complex creation failed." << std::endl; + exit(-1); + } + } + return stree; +} + int main(int argc, char **argv) { + std::string weight_file; std::string off_file_points; std::string output_file_diag; bool exact_version = false; @@ -38,48 +101,34 @@ int main(int argc, char **argv) { int coeff_field_characteristic; Filtration_value min_persistence; - program_options(argc, argv, off_file_points, exact_version, fast_version, output_file_diag, alpha_square_max_value, - coeff_field_characteristic, min_persistence); + program_options(argc, argv, off_file_points, exact_version, fast_version, weight_file, output_file_diag, + alpha_square_max_value, coeff_field_characteristic, min_persistence); if ((exact_version) && (fast_version)) { std::cerr << "You cannot set the exact and the fast version." << std::endl; exit(-1); } - Simplex_tree simplex; + Simplex_tree stree; if (fast_version) { // WARNING : CGAL::Epick_d is fast but not safe (unlike CGAL::Epeck_d) // (i.e. when the points are on a grid) using Fast_kernel = CGAL::Epick_d; - - // Init of an alpha complex from an OFF file - Gudhi::alpha_complex::Alpha_complex alpha_complex_from_file(off_file_points); - - if (!alpha_complex_from_file.create_complex(simplex, alpha_square_max_value)) { - std::cerr << "Fast Alpha complex simplicial complex creation failed." << std::endl; - exit(-1); - } + stree = create_simplex_tree(off_file_points, weight_file, exact_version, alpha_square_max_value); } else { using Kernel = CGAL::Epeck_d; - - // Init of an alpha complex from an OFF file - Gudhi::alpha_complex::Alpha_complex alpha_complex_from_file(off_file_points); - - if (!alpha_complex_from_file.create_complex(simplex, alpha_square_max_value, exact_version)) { - std::cerr << "Alpha complex simplicial complex creation failed." << std::endl; - exit(-1); - } + stree = create_simplex_tree(off_file_points, weight_file, exact_version, alpha_square_max_value); } // ---------------------------------------------------------------------------- // Display information about the alpha complex // ---------------------------------------------------------------------------- - std::clog << "Simplicial complex is of dimension " << simplex.dimension() << " - " << simplex.num_simplices() - << " simplices - " << simplex.num_vertices() << " vertices." << std::endl; + std::clog << "Simplicial complex is of dimension " << stree.dimension() << " - " << stree.num_simplices() + << " simplices - " << stree.num_vertices() << " vertices." << std::endl; - std::clog << "Simplex_tree dim: " << simplex.dimension() << std::endl; + std::clog << "Simplex_tree dim: " << stree.dimension() << std::endl; // Compute the persistence diagram of the complex Gudhi::persistent_cohomology::Persistent_cohomology pcoh( - simplex); + stree); // initializes the coefficient field for homology pcoh.init_coefficients(coeff_field_characteristic); @@ -98,7 +147,7 @@ int main(int argc, char **argv) { } void program_options(int argc, char *argv[], std::string &off_file_points, bool &exact, bool &fast, - std::string &output_file_diag, Filtration_value &alpha_square_max_value, + std::string &weight_file, std::string &output_file_diag, Filtration_value &alpha_square_max_value, int &coeff_field_characteristic, Filtration_value &min_persistence) { namespace po = boost::program_options; po::options_description hidden("Hidden options"); @@ -111,6 +160,8 @@ void program_options(int argc, char *argv[], std::string &off_file_points, bool "To activate exact version of Alpha complex (default is false, not available if fast is set)")( "fast,f", po::bool_switch(&fast), "To activate fast version of Alpha complex (default is false, not available if exact is set)")( + "weight-file,w", po::value(&weight_file)->default_value(std::string()), + "Name of file containing a point weights. Format is one weight per line:\n W1\n ...\n Wn ")( "output-file,o", po::value(&output_file_diag)->default_value(std::string()), "Name of file in which the persistence diagram is written. Default print in std::clog")( "max-alpha-square-value,r", po::value(&alpha_square_max_value) -- cgit v1.2.3 From c1b317640f81fefff7f7a774b239957fb632fa43 Mon Sep 17 00:00:00 2001 From: ROUVREAU Vincent Date: Wed, 30 Sep 2020 09:18:36 +0200 Subject: code review: emplace_back versus push_back --- .../Weighted_alpha_complex_3d_from_points.cpp | 10 +++++----- .../example/Weighted_alpha_complex_from_points.cpp | 10 +++++----- src/Alpha_complex/include/gudhi/Alpha_complex.h | 6 +++--- src/Alpha_complex/test/Alpha_kernel_d_unit_test.cpp | 20 ++++++++++---------- .../utilities/alpha_complex_persistence.cpp | 2 +- 5 files changed, 24 insertions(+), 24 deletions(-) (limited to 'src/Alpha_complex/utilities/alpha_complex_persistence.cpp') diff --git a/src/Alpha_complex/example/Weighted_alpha_complex_3d_from_points.cpp b/src/Alpha_complex/example/Weighted_alpha_complex_3d_from_points.cpp index a8d00272..09431d6e 100644 --- a/src/Alpha_complex/example/Weighted_alpha_complex_3d_from_points.cpp +++ b/src/Alpha_complex/example/Weighted_alpha_complex_3d_from_points.cpp @@ -18,11 +18,11 @@ int main(int argc, char **argv) { // Init of a list of points and weights from a small molecule // ---------------------------------------------------------------------------- std::vector weighted_points; - weighted_points.push_back(Weighted_point(Bare_point(1, -1, -1), 4.)); - weighted_points.push_back(Weighted_point(Bare_point(-1, 1, -1), 4.)); - weighted_points.push_back(Weighted_point(Bare_point(-1, -1, 1), 4.)); - weighted_points.push_back(Weighted_point(Bare_point(1, 1, 1), 4.)); - weighted_points.push_back(Weighted_point(Bare_point(2, 2, 2), 1.)); + weighted_points.emplace_back(Bare_point(1, -1, -1), 4.); + weighted_points.emplace_back(Bare_point(-1, 1, -1), 4.); + weighted_points.emplace_back(Bare_point(-1, -1, 1), 4.); + weighted_points.emplace_back(Bare_point(1, 1, 1), 4.); + weighted_points.emplace_back(Bare_point(2, 2, 2), 1.); // ---------------------------------------------------------------------------- // Init of an alpha complex from the list of points diff --git a/src/Alpha_complex/example/Weighted_alpha_complex_from_points.cpp b/src/Alpha_complex/example/Weighted_alpha_complex_from_points.cpp index d49d3e93..d1f3e436 100644 --- a/src/Alpha_complex/example/Weighted_alpha_complex_from_points.cpp +++ b/src/Alpha_complex/example/Weighted_alpha_complex_from_points.cpp @@ -18,11 +18,11 @@ int main() { // Init of a list of points and weights from a small molecule // ---------------------------------------------------------------------------- Vector_of_points points; - points.push_back(Weighted_point(Bare_point(1, -1, -1), 4.)); - points.push_back(Weighted_point(Bare_point(-1, 1, -1), 4.)); - points.push_back(Weighted_point(Bare_point(-1, -1, 1), 4.)); - points.push_back(Weighted_point(Bare_point(1, 1, 1), 4.)); - points.push_back(Weighted_point(Bare_point(2, 2, 2), 1.)); + points.emplace_back(Bare_point(1, -1, -1), 4.); + points.emplace_back(Bare_point(-1, 1, -1), 4.); + points.emplace_back(Bare_point(-1, -1, 1), 4.); + points.emplace_back(Bare_point(1, 1, 1), 4.); + points.emplace_back(Bare_point(2, 2, 2), 1.); // ---------------------------------------------------------------------------- // Init of an alpha complex from the list of points diff --git a/src/Alpha_complex/include/gudhi/Alpha_complex.h b/src/Alpha_complex/include/gudhi/Alpha_complex.h index a692dbc6..4f8e3d54 100644 --- a/src/Alpha_complex/include/gudhi/Alpha_complex.h +++ b/src/Alpha_complex/include/gudhi/Alpha_complex.h @@ -310,7 +310,7 @@ class Alpha_complex { thread_local std::vector v; v.clear(); for (auto vertex : cplx.simplex_vertex_range(s)) - v.push_back(get_point_(vertex)); + v.emplace_back(get_point_(vertex)); cache_.emplace_back(kernel_.get_sphere(v.cbegin(), v.cend())); } return cache_[k]; @@ -326,7 +326,7 @@ class Alpha_complex { thread_local std::vector v; v.clear(); for (auto vertex : cplx.simplex_vertex_range(s)) - v.push_back(get_point_(vertex)); + v.emplace_back(get_point_(vertex)); return kernel_.get_squared_radius(v.cbegin(), v.cend()); } @@ -394,7 +394,7 @@ class Alpha_complex { std::clog << " " << (*vit)->data(); #endif // DEBUG_TRACES // Vector of vertex construction for simplex_tree structure - vertexVector.push_back((*vit)->data()); + vertexVector.emplace_back((*vit)->data()); } } #ifdef DEBUG_TRACES diff --git a/src/Alpha_complex/test/Alpha_kernel_d_unit_test.cpp b/src/Alpha_complex/test/Alpha_kernel_d_unit_test.cpp index 192834b3..6eae103d 100644 --- a/src/Alpha_complex/test/Alpha_kernel_d_unit_test.cpp +++ b/src/Alpha_complex/test/Alpha_kernel_d_unit_test.cpp @@ -63,11 +63,11 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(Alpha_kernel_d_sphere, TestedKernel, list_of_kerne using Point_d = typename Unweighted_kernel::Point_d; std::vector unw_pts; - unw_pts.push_back(Point_d(p0.begin(), p0.end())); - unw_pts.push_back(Point_d(p1.begin(), p1.end())); - unw_pts.push_back(Point_d(p2.begin(), p2.end())); - unw_pts.push_back(Point_d(p3.begin(), p3.end())); - unw_pts.push_back(Point_d(p4.begin(), p4.end())); + unw_pts.emplace_back(p0.begin(), p0.end()); + unw_pts.emplace_back(p1.begin(), p1.end()); + unw_pts.emplace_back(p2.begin(), p2.end()); + unw_pts.emplace_back(p3.begin(), p3.end()); + unw_pts.emplace_back(p4.begin(), p4.end()); Unweighted_kernel kernel; auto unw_sphere = kernel.get_sphere(unw_pts.cbegin(), unw_pts.cend()); @@ -79,11 +79,11 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(Alpha_kernel_d_sphere, TestedKernel, list_of_kerne using Weighted_point_d = typename Weighted_kernel::Weighted_point_d; using Bare_point_d = typename Weighted_kernel::Bare_point_d; std::vector w_pts; - w_pts.push_back(Weighted_point_d(Bare_point_d(p0.begin(), p0.end()), 0.)); - w_pts.push_back(Weighted_point_d(Bare_point_d(p1.begin(), p1.end()), 0.)); - w_pts.push_back(Weighted_point_d(Bare_point_d(p2.begin(), p2.end()), 0.)); - w_pts.push_back(Weighted_point_d(Bare_point_d(p3.begin(), p3.end()), 0.)); - w_pts.push_back(Weighted_point_d(Bare_point_d(p4.begin(), p4.end()), 0.)); + w_pts.emplace_back(Bare_point_d(p0.begin(), p0.end()), 0.); + w_pts.emplace_back(Bare_point_d(p1.begin(), p1.end()), 0.); + w_pts.emplace_back(Bare_point_d(p2.begin(), p2.end()), 0.); + w_pts.emplace_back(Bare_point_d(p3.begin(), p3.end()), 0.); + w_pts.emplace_back(Bare_point_d(p4.begin(), p4.end()), 0.); Weighted_kernel w_kernel; auto w_sphere = w_kernel.get_sphere(w_pts.cbegin(), w_pts.cend()); diff --git a/src/Alpha_complex/utilities/alpha_complex_persistence.cpp b/src/Alpha_complex/utilities/alpha_complex_persistence.cpp index e86b34e2..3ce7b440 100644 --- a/src/Alpha_complex/utilities/alpha_complex_persistence.cpp +++ b/src/Alpha_complex/utilities/alpha_complex_persistence.cpp @@ -50,7 +50,7 @@ std::vector read_weight_file(const std::string &weight_file) { double weight = 0.0; // Attempt read the weight in a double format, return false if it fails while (weights_ifstr >> weight) { - weights.push_back(weight); + weights.emplace_back(weight); } } else { std::cerr << "Unable to read weights file " << weight_file << std::endl; -- cgit v1.2.3 From 74143b372d3d93fc2139974fa85906047acee8a7 Mon Sep 17 00:00:00 2001 From: ROUVREAU Vincent Date: Mon, 5 Oct 2020 10:33:23 +0200 Subject: code review: push_back versus emplace_back --- src/Alpha_complex/include/gudhi/Alpha_complex.h | 6 +++--- src/Alpha_complex/utilities/alpha_complex_persistence.cpp | 2 +- 2 files changed, 4 insertions(+), 4 deletions(-) (limited to 'src/Alpha_complex/utilities/alpha_complex_persistence.cpp') diff --git a/src/Alpha_complex/include/gudhi/Alpha_complex.h b/src/Alpha_complex/include/gudhi/Alpha_complex.h index 4f8e3d54..a692dbc6 100644 --- a/src/Alpha_complex/include/gudhi/Alpha_complex.h +++ b/src/Alpha_complex/include/gudhi/Alpha_complex.h @@ -310,7 +310,7 @@ class Alpha_complex { thread_local std::vector v; v.clear(); for (auto vertex : cplx.simplex_vertex_range(s)) - v.emplace_back(get_point_(vertex)); + v.push_back(get_point_(vertex)); cache_.emplace_back(kernel_.get_sphere(v.cbegin(), v.cend())); } return cache_[k]; @@ -326,7 +326,7 @@ class Alpha_complex { thread_local std::vector v; v.clear(); for (auto vertex : cplx.simplex_vertex_range(s)) - v.emplace_back(get_point_(vertex)); + v.push_back(get_point_(vertex)); return kernel_.get_squared_radius(v.cbegin(), v.cend()); } @@ -394,7 +394,7 @@ class Alpha_complex { std::clog << " " << (*vit)->data(); #endif // DEBUG_TRACES // Vector of vertex construction for simplex_tree structure - vertexVector.emplace_back((*vit)->data()); + vertexVector.push_back((*vit)->data()); } } #ifdef DEBUG_TRACES diff --git a/src/Alpha_complex/utilities/alpha_complex_persistence.cpp b/src/Alpha_complex/utilities/alpha_complex_persistence.cpp index 3ce7b440..e86b34e2 100644 --- a/src/Alpha_complex/utilities/alpha_complex_persistence.cpp +++ b/src/Alpha_complex/utilities/alpha_complex_persistence.cpp @@ -50,7 +50,7 @@ std::vector read_weight_file(const std::string &weight_file) { double weight = 0.0; // Attempt read the weight in a double format, return false if it fails while (weights_ifstr >> weight) { - weights.emplace_back(weight); + weights.push_back(weight); } } else { std::cerr << "Unable to read weights file " << weight_file << std::endl; -- cgit v1.2.3