diff options
Diffstat (limited to 'src/Alpha_complex')
21 files changed, 399 insertions, 286 deletions
diff --git a/src/Alpha_complex/benchmark/Alpha_complex_3d_benchmark.cpp b/src/Alpha_complex/benchmark/Alpha_complex_3d_benchmark.cpp index 99ad94b9..e7d85686 100644 --- a/src/Alpha_complex/benchmark/Alpha_complex_3d_benchmark.cpp +++ b/src/Alpha_complex/benchmark/Alpha_complex_3d_benchmark.cpp @@ -19,7 +19,7 @@ std::ofstream results_csv("results.csv"); template <typename Kernel> void benchmark_points_on_torus_dD(const std::string& msg) { - std::cout << "+ " << msg << std::endl; + std::clog << "+ " << msg << std::endl; results_csv << "\"" << msg << "\";" << std::endl; results_csv << "\"nb_points\";" @@ -29,7 +29,7 @@ void benchmark_points_on_torus_dD(const std::string& msg) { using K = CGAL::Epick_d<CGAL::Dimension_tag<3>>; for (int nb_points = 1000; nb_points <= 125000; nb_points *= 5) { - std::cout << " Alpha complex dD on torus with " << nb_points << " points." << std::endl; + std::clog << " Alpha complex dD on torus with " << nb_points << " points." << std::endl; std::vector<K::Point_d> points_on_torus = Gudhi::generate_points_on_torus_3D<K>(nb_points, 1.0, 0.5); std::vector<typename Kernel::Point_d> points; @@ -41,26 +41,26 @@ void benchmark_points_on_torus_dD(const std::string& msg) { ac_create_clock.begin(); Gudhi::alpha_complex::Alpha_complex<Kernel> alpha_complex_from_points(points); ac_create_clock.end(); - std::cout << ac_create_clock; + std::clog << ac_create_clock; Gudhi::Simplex_tree<> complex; Gudhi::Clock st_create_clock(" benchmark_points_on_torus_dD - complex creation"); st_create_clock.begin(); alpha_complex_from_points.create_complex(complex); st_create_clock.end(); - std::cout << st_create_clock; + std::clog << st_create_clock; results_csv << nb_points << ";" << complex.num_simplices() << ";" << ac_create_clock.num_seconds() << ";" << st_create_clock.num_seconds() << ";" << std::endl; - std::cout << " benchmark_points_on_torus_dD - nb simplices = " << complex.num_simplices() << std::endl; + std::clog << " benchmark_points_on_torus_dD - nb simplices = " << complex.num_simplices() << std::endl; } } template <typename Alpha_complex_3d> void benchmark_points_on_torus_3D(const std::string& msg) { using K = CGAL::Epick_d<CGAL::Dimension_tag<3>>; - std::cout << "+ " << msg << std::endl; + std::clog << "+ " << msg << std::endl; results_csv << "\"" << msg << "\";" << std::endl; results_csv << "\"nb_points\";" @@ -69,7 +69,7 @@ void benchmark_points_on_torus_3D(const std::string& msg) { << "\"complex_creation_time(sec.)\";" << std::endl; for (int nb_points = 1000; nb_points <= 125000; nb_points *= 5) { - std::cout << " Alpha complex 3d on torus with " << nb_points << " points." << std::endl; + std::clog << " Alpha complex 3d on torus with " << nb_points << " points." << std::endl; std::vector<K::Point_d> points_on_torus = Gudhi::generate_points_on_torus_3D<K>(nb_points, 1.0, 0.5); std::vector<typename Alpha_complex_3d::Point_3> points; @@ -81,19 +81,19 @@ void benchmark_points_on_torus_3D(const std::string& msg) { ac_create_clock.begin(); Alpha_complex_3d alpha_complex_from_points(points); ac_create_clock.end(); - std::cout << ac_create_clock; + std::clog << ac_create_clock; Gudhi::Simplex_tree<> complex; Gudhi::Clock st_create_clock(" benchmark_points_on_torus_3D - complex creation"); st_create_clock.begin(); alpha_complex_from_points.create_complex(complex); st_create_clock.end(); - std::cout << st_create_clock; + std::clog << st_create_clock; results_csv << nb_points << ";" << complex.num_simplices() << ";" << ac_create_clock.num_seconds() << ";" << st_create_clock.num_seconds() << ";" << std::endl; - std::cout << " benchmark_points_on_torus_3D - nb simplices = " << complex.num_simplices() << std::endl; + std::clog << " benchmark_points_on_torus_3D - nb simplices = " << complex.num_simplices() << std::endl; } } @@ -101,7 +101,7 @@ template <typename Weighted_alpha_complex_3d> void benchmark_weighted_points_on_torus_3D(const std::string& msg) { using K = CGAL::Epick_d<CGAL::Dimension_tag<3>>; - std::cout << "+ " << msg << std::endl; + std::clog << "+ " << msg << std::endl; results_csv << "\"" << msg << "\";" << std::endl; results_csv << "\"nb_points\";" @@ -112,7 +112,7 @@ void benchmark_weighted_points_on_torus_3D(const std::string& msg) { CGAL::Random random(8); for (int nb_points = 1000; nb_points <= 125000; nb_points *= 5) { - std::cout << " Alpha complex 3d on torus with " << nb_points << " points." << std::endl; + std::clog << " Alpha complex 3d on torus with " << nb_points << " points." << std::endl; std::vector<K::Point_d> points_on_torus = Gudhi::generate_points_on_torus_3D<K>(nb_points, 1.0, 0.5); using Point = typename Weighted_alpha_complex_3d::Bare_point_3; @@ -128,25 +128,25 @@ void benchmark_weighted_points_on_torus_3D(const std::string& msg) { ac_create_clock.begin(); Weighted_alpha_complex_3d alpha_complex_from_points(points); ac_create_clock.end(); - std::cout << ac_create_clock; + std::clog << ac_create_clock; Gudhi::Simplex_tree<> complex; Gudhi::Clock st_create_clock(" benchmark_weighted_points_on_torus_3D - complex creation"); st_create_clock.begin(); alpha_complex_from_points.create_complex(complex); st_create_clock.end(); - std::cout << st_create_clock; + std::clog << st_create_clock; results_csv << nb_points << ";" << complex.num_simplices() << ";" << ac_create_clock.num_seconds() << ";" << st_create_clock.num_seconds() << ";" << std::endl; - std::cout << " benchmark_weighted_points_on_torus_3D - nb simplices = " << complex.num_simplices() << std::endl; + std::clog << " benchmark_weighted_points_on_torus_3D - nb simplices = " << complex.num_simplices() << std::endl; } } template <typename Periodic_alpha_complex_3d> void benchmark_periodic_points(const std::string& msg) { - std::cout << "+ " << msg << std::endl; + std::clog << "+ " << msg << std::endl; results_csv << "\"" << msg << "\";" << std::endl; results_csv << "\"nb_points\";" @@ -157,7 +157,7 @@ void benchmark_periodic_points(const std::string& msg) { CGAL::Random random(8); for (double nb_points = 10.; nb_points <= 40.; nb_points += 10.) { - std::cout << " Periodic alpha complex 3d with " << nb_points * nb_points * nb_points << " points." << std::endl; + std::clog << " Periodic alpha complex 3d with " << nb_points * nb_points * nb_points << " points." << std::endl; using Point = typename Periodic_alpha_complex_3d::Point_3; std::vector<Point> points; @@ -174,25 +174,25 @@ void benchmark_periodic_points(const std::string& msg) { ac_create_clock.begin(); Periodic_alpha_complex_3d alpha_complex_from_points(points, 0., 0., 0., nb_points, nb_points, nb_points); ac_create_clock.end(); - std::cout << ac_create_clock; + std::clog << ac_create_clock; Gudhi::Simplex_tree<> complex; Gudhi::Clock st_create_clock(" benchmark_periodic_points - complex creation"); st_create_clock.begin(); alpha_complex_from_points.create_complex(complex); st_create_clock.end(); - std::cout << st_create_clock; + std::clog << st_create_clock; results_csv << nb_points * nb_points * nb_points << ";" << complex.num_simplices() << ";" << ac_create_clock.num_seconds() << ";" << st_create_clock.num_seconds() << ";" << std::endl; - std::cout << " benchmark_periodic_points - nb simplices = " << complex.num_simplices() << std::endl; + std::clog << " benchmark_periodic_points - nb simplices = " << complex.num_simplices() << std::endl; } } template <typename Weighted_periodic_alpha_complex_3d> void benchmark_weighted_periodic_points(const std::string& msg) { - std::cout << "+ " << msg << std::endl; + std::clog << "+ " << msg << std::endl; results_csv << "\"" << msg << "\";" << std::endl; results_csv << "\"nb_points\";" @@ -203,7 +203,7 @@ void benchmark_weighted_periodic_points(const std::string& msg) { CGAL::Random random(8); for (double nb_points = 10.; nb_points <= 40.; nb_points += 10.) { - std::cout << " Weighted periodic alpha complex 3d with " << nb_points * nb_points * nb_points << " points." + std::clog << " Weighted periodic alpha complex 3d with " << nb_points * nb_points * nb_points << " points." << std::endl; using Point = typename Weighted_periodic_alpha_complex_3d::Bare_point_3; @@ -224,19 +224,19 @@ void benchmark_weighted_periodic_points(const std::string& msg) { ac_create_clock.begin(); Weighted_periodic_alpha_complex_3d alpha_complex_from_points(points, 0., 0., 0., nb_points, nb_points, nb_points); ac_create_clock.end(); - std::cout << ac_create_clock; + std::clog << ac_create_clock; Gudhi::Simplex_tree<> complex; Gudhi::Clock st_create_clock(" benchmark_weighted_periodic_points - complex creation"); st_create_clock.begin(); alpha_complex_from_points.create_complex(complex); st_create_clock.end(); - std::cout << st_create_clock; + std::clog << st_create_clock; results_csv << nb_points * nb_points * nb_points << ";" << complex.num_simplices() << ";" << ac_create_clock.num_seconds() << ";" << st_create_clock.num_seconds() << ";" << std::endl; - std::cout << " benchmark_weighted_periodic_points - nb simplices = " << complex.num_simplices() << std::endl; + std::clog << " benchmark_weighted_periodic_points - nb simplices = " << complex.num_simplices() << std::endl; } } diff --git a/src/Alpha_complex/concept/SimplicialComplexForAlpha.h b/src/Alpha_complex/concept/SimplicialComplexForAlpha.h index 1c6c3b0c..c20c3201 100644 --- a/src/Alpha_complex/concept/SimplicialComplexForAlpha.h +++ b/src/Alpha_complex/concept/SimplicialComplexForAlpha.h @@ -72,6 +72,24 @@ struct SimplicialComplexForAlpha { /** \brief Return type of an insertion of a simplex */ typedef unspecified Insertion_result_type; + + /** \name Map interface + * Conceptually a `std::unordered_map<Simplex_handle,std::size_t>`. + * @{ */ + /** \brief Data stored for each simplex. + * + * Must be an integer type. */ + typedef unspecified Simplex_key; + /** \brief Returns a constant dummy number that is either negative, + * or at least as large as the number of simplices. Suggested value: -1. */ + Simplex_key null_key (); + /** \brief Returns the number stored for a simplex by `assign_key()`. + * + * If `assign_key()` has not been called, it must return `null_key()`. */ + Simplex_key key ( Simplex_handle sh ); + /** \brief Store a number for a simplex, which can later be retrieved with `key()`. */ + void assign_key(Simplex_handle sh, Simplex_key n); + /** @} */ }; } // namespace alpha_complex diff --git a/src/Alpha_complex/doc/Intro_alpha_complex.h b/src/Alpha_complex/doc/Intro_alpha_complex.h index a8b1a106..60da7169 100644 --- a/src/Alpha_complex/doc/Intro_alpha_complex.h +++ b/src/Alpha_complex/doc/Intro_alpha_complex.h @@ -46,16 +46,17 @@ namespace alpha_complex { * \cite cgal:s-gkd-19b from CGAL as template parameter. * * \remark - * - When an \f$\alpha\f$-complex is constructed with an infinite value of \f$ \alpha^2 \f$, the complex is a Delaunay - * complex (with special filtration values). + * - When the simplicial complex is constructed with an infinite value of \f$ \alpha^2 \f$, the complex is a Delaunay + * complex with special filtration values. The Delaunay complex without filtration values is also available by passing + * `default_filtration_value=true` to `Alpha_complex::create_complex`. * - For people only interested in the topology of the \ref alpha_complex (for instance persistence), * \ref alpha_complex is equivalent to the \ref cech_complex and much smaller if you do not bound the radii. * \ref cech_complex can still make sense in higher dimension precisely because you can bound the radii. - * - Using the default `CGAL::Epeck_d` makes the construction safe. If you pass exact=true to create_complex, the - * filtration values are the exact ones converted to the filtration value type of the simplicial complex. This can be - * very slow. If you pass exact=false (the default), the filtration values are only guaranteed to have a small - * multiplicative error compared to the exact value, see <code><a class="el" target="_blank" - * href="https://doc.cgal.org/latest/Number_types/classCGAL_1_1Lazy__exact__nt.html"> + * - Using the default `CGAL::Epeck_d` makes the construction safe. If you pass `exact=true` to + * `Alpha_complex::create_complex`, the filtration values are the exact ones converted to the filtration value type of + * the simplicial complex. This can be very slow. If you pass `exact=false` (the default), the filtration values are + * only guaranteed to have a small multiplicative error compared to the exact value, see <code> + * <a class="el" target="_blank" href="https://doc.cgal.org/latest/Number_types/classCGAL_1_1Lazy__exact__nt.html"> * CGAL::Lazy_exact_nt<NT>::set_relative_precision_of_to_double</a></code> for details. A drawback, when computing * persistence, is that an empty exact interval [10^12,10^12] may become a non-empty approximate interval * [10^12,10^12+10^6]. Using `CGAL::Epick_d` makes the computations slightly faster, and the combinatorics are still diff --git a/src/Alpha_complex/example/Alpha_complex_3d_from_points.cpp b/src/Alpha_complex/example/Alpha_complex_3d_from_points.cpp index 0e359a27..a2c85138 100644 --- a/src/Alpha_complex/example/Alpha_complex_3d_from_points.cpp +++ b/src/Alpha_complex/example/Alpha_complex_3d_from_points.cpp @@ -38,18 +38,18 @@ int main(int argc, char **argv) { // ---------------------------------------------------------------------------- // Display information about the alpha complex // ---------------------------------------------------------------------------- - std::cout << "Alpha complex is of dimension " << simplex.dimension() << " - " << simplex.num_simplices() + std::clog << "Alpha complex is of dimension " << simplex.dimension() << " - " << simplex.num_simplices() << " simplices - " << simplex.num_vertices() << " vertices." << std::endl; - std::cout << "Iterator on alpha complex simplices in the filtration order, with [filtration value]:" << std::endl; + std::clog << "Iterator on alpha complex simplices in the filtration order, with [filtration value]:" << std::endl; for (auto f_simplex : simplex.filtration_simplex_range()) { - std::cout << " ( "; + std::clog << " ( "; for (auto vertex : simplex.simplex_vertex_range(f_simplex)) { - std::cout << vertex << " "; + std::clog << vertex << " "; } - std::cout << ") -> " + std::clog << ") -> " << "[" << simplex.filtration(f_simplex) << "] "; - std::cout << std::endl; + std::clog << std::endl; } } return 0; diff --git a/src/Alpha_complex/example/Alpha_complex_from_off.cpp b/src/Alpha_complex/example/Alpha_complex_from_off.cpp index 220a66de..dba1710e 100644 --- a/src/Alpha_complex/example/Alpha_complex_from_off.cpp +++ b/src/Alpha_complex/example/Alpha_complex_from_off.cpp @@ -30,7 +30,7 @@ int main(int argc, char **argv) { ouput_file_stream.open(std::string(argv[3])); streambuffer = ouput_file_stream.rdbuf(); } else { - streambuffer = std::cout.rdbuf(); + streambuffer = std::clog.rdbuf(); } Gudhi::Simplex_tree<> simplex; diff --git a/src/Alpha_complex/example/Alpha_complex_from_points.cpp b/src/Alpha_complex/example/Alpha_complex_from_points.cpp index 6526ca3a..c79535bf 100644 --- a/src/Alpha_complex/example/Alpha_complex_from_points.cpp +++ b/src/Alpha_complex/example/Alpha_complex_from_points.cpp @@ -35,18 +35,18 @@ int main() { // ---------------------------------------------------------------------------- // Display information about the alpha complex // ---------------------------------------------------------------------------- - std::cout << "Alpha complex is of dimension " << simplex.dimension() << + std::clog << "Alpha complex is of dimension " << simplex.dimension() << " - " << simplex.num_simplices() << " simplices - " << simplex.num_vertices() << " vertices." << std::endl; - std::cout << "Iterator on alpha complex simplices in the filtration order, with [filtration value]:" << std::endl; + std::clog << "Iterator on alpha complex simplices in the filtration order, with [filtration value]:" << std::endl; for (auto f_simplex : simplex.filtration_simplex_range()) { - std::cout << " ( "; + std::clog << " ( "; for (auto vertex : simplex.simplex_vertex_range(f_simplex)) { - std::cout << vertex << " "; + std::clog << vertex << " "; } - std::cout << ") -> " << "[" << simplex.filtration(f_simplex) << "] "; - std::cout << std::endl; + std::clog << ") -> " << "[" << simplex.filtration(f_simplex) << "] "; + std::clog << std::endl; } } return 0; diff --git a/src/Alpha_complex/example/CMakeLists.txt b/src/Alpha_complex/example/CMakeLists.txt index b0337934..2eecd50c 100644 --- a/src/Alpha_complex/example/CMakeLists.txt +++ b/src/Alpha_complex/example/CMakeLists.txt @@ -32,14 +32,18 @@ if (DIFF_PATH) add_test(Alpha_complex_example_from_off_60_diff_files ${DIFF_PATH} ${CMAKE_CURRENT_BINARY_DIR}/alphaoffreader_result_60.txt ${CMAKE_CURRENT_BINARY_DIR}/alphaoffreader_for_doc_60.txt) + set_tests_properties(Alpha_complex_example_from_off_60_diff_files PROPERTIES DEPENDS Alpha_complex_example_from_off_60) add_test(Alpha_complex_example_from_off_32_diff_files ${DIFF_PATH} ${CMAKE_CURRENT_BINARY_DIR}/alphaoffreader_result_32.txt ${CMAKE_CURRENT_BINARY_DIR}/alphaoffreader_for_doc_32.txt) + set_tests_properties(Alpha_complex_example_from_off_32_diff_files PROPERTIES DEPENDS Alpha_complex_example_from_off_32) add_test(Alpha_complex_example_fast_from_off_60_diff_files ${DIFF_PATH} ${CMAKE_CURRENT_BINARY_DIR}/fastalphaoffreader_result_60.txt ${CMAKE_CURRENT_BINARY_DIR}/alphaoffreader_for_doc_60.txt) + set_tests_properties(Alpha_complex_example_fast_from_off_60_diff_files PROPERTIES DEPENDS Alpha_complex_example_fast_from_off_60) add_test(Alpha_complex_example_fast_from_off_32_diff_files ${DIFF_PATH} ${CMAKE_CURRENT_BINARY_DIR}/fastalphaoffreader_result_32.txt ${CMAKE_CURRENT_BINARY_DIR}/alphaoffreader_for_doc_32.txt) - endif() + set_tests_properties(Alpha_complex_example_fast_from_off_32_diff_files PROPERTIES DEPENDS Alpha_complex_example_fast_from_off_32) +endif() add_executable ( Alpha_complex_example_weighted_3d_from_points Weighted_alpha_complex_3d_from_points.cpp ) target_link_libraries(Alpha_complex_example_weighted_3d_from_points ${CGAL_LIBRARY}) diff --git a/src/Alpha_complex/example/Fast_alpha_complex_from_off.cpp b/src/Alpha_complex/example/Fast_alpha_complex_from_off.cpp index f181005a..64728470 100644 --- a/src/Alpha_complex/example/Fast_alpha_complex_from_off.cpp +++ b/src/Alpha_complex/example/Fast_alpha_complex_from_off.cpp @@ -35,7 +35,7 @@ int main(int argc, char **argv) { ouput_file_stream.open(std::string(argv[3])); streambuffer = ouput_file_stream.rdbuf(); } else { - streambuffer = std::cout.rdbuf(); + streambuffer = std::clog.rdbuf(); } Gudhi::Simplex_tree<> simplex; 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 fcf80802..c044194e 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 @@ -34,18 +34,18 @@ int main(int argc, char **argv) { // ---------------------------------------------------------------------------- // Display information about the alpha complex // ---------------------------------------------------------------------------- - std::cout << "Alpha complex is of dimension " << simplex.dimension() << " - " << simplex.num_simplices() + std::clog << "Alpha complex is of dimension " << simplex.dimension() << " - " << simplex.num_simplices() << " simplices - " << simplex.num_vertices() << " vertices." << std::endl; - std::cout << "Iterator on alpha complex simplices in the filtration order, with [filtration value]:" << std::endl; + std::clog << "Iterator on alpha complex simplices in the filtration order, with [filtration value]:" << std::endl; for (auto f_simplex : simplex.filtration_simplex_range()) { - std::cout << " ( "; + std::clog << " ( "; for (auto vertex : simplex.simplex_vertex_range(f_simplex)) { - std::cout << vertex << " "; + std::clog << vertex << " "; } - std::cout << ") -> " + std::clog << ") -> " << "[" << simplex.filtration(f_simplex) << "] "; - std::cout << std::endl; + std::clog << std::endl; } } return 0; diff --git a/src/Alpha_complex/include/gudhi/Alpha_complex.h b/src/Alpha_complex/include/gudhi/Alpha_complex.h index f2a05e95..ba91998d 100644 --- a/src/Alpha_complex/include/gudhi/Alpha_complex.h +++ b/src/Alpha_complex/include/gudhi/Alpha_complex.h @@ -112,9 +112,6 @@ class Alpha_complex { typedef typename Kernel::Side_of_bounded_sphere_d Is_Gabriel; typedef typename Kernel::Point_dimension_d Point_Dimension; - // Type required to compute squared radius, or side of bounded sphere on a vector of points. - typedef typename std::vector<Point_d> Vector_of_CGAL_points; - // Vertex_iterator type from CGAL. typedef typename Delaunay_triangulation::Vertex_iterator CGAL_vertex_iterator; @@ -124,6 +121,9 @@ class Alpha_complex { // Structure to switch from simplex tree vertex handle to CGAL vertex iterator. typedef typename std::vector< CGAL_vertex_iterator > Vector_vertex_iterator; + // Numeric type of coordinates in the kernel + typedef typename Kernel::FT FT; + private: /** \brief Vertex iterator vector to switch from simplex tree vertex handle to CGAL vertex iterator. * Vertex handles are inserted sequentially, starting at 0.*/ @@ -132,6 +132,8 @@ class Alpha_complex { Delaunay_triangulation* triangulation_; /** \brief Kernel for triangulation_ functions access.*/ Kernel kernel_; + /** \brief Cache for geometric constructions: circumcenter and squared radius of a simplex.*/ + std::vector<std::pair<Point_d, FT>> cache_, old_cache_; public: /** \brief Alpha_complex constructor from an OFF file name. @@ -237,7 +239,7 @@ class Alpha_complex { for (CGAL_vertex_iterator vit = triangulation_->vertices_begin(); vit != triangulation_->vertices_end(); ++vit) { if (!triangulation_->is_infinite(*vit)) { #ifdef DEBUG_TRACES - std::cout << "Vertex insertion - " << vit->data() << " -> " << vit->point() << std::endl; + std::clog << "Vertex insertion - " << vit->data() << " -> " << vit->point() << std::endl; #endif // DEBUG_TRACES vertex_handle_to_iterator_[vit->data()] = vit; } @@ -246,18 +248,66 @@ class Alpha_complex { } } + /** \brief get_point_ returns the point corresponding to the vertex given as parameter. + * Only for internal use for faster access. + * + * @param[in] vertex Vertex handle of the point to retrieve. + * @return The point found. + */ + const Point_d& get_point_(std::size_t vertex) const { + return vertex_handle_to_iterator_[vertex]->point(); + } + + /// Return a reference to the circumcenter and circumradius, writing them in the cache if necessary. + template<class SimplicialComplexForAlpha> + auto& get_cache(SimplicialComplexForAlpha& cplx, typename SimplicialComplexForAlpha::Simplex_handle s) { + auto k = cplx.key(s); + if(k==cplx.null_key()){ + k = cache_.size(); + cplx.assign_key(s, k); + // Using a transform_range is slower, currently. + thread_local std::vector<Point_d> v; + v.clear(); + for (auto vertex : cplx.simplex_vertex_range(s)) + v.push_back(get_point_(vertex)); + Point_d c = kernel_.construct_circumcenter_d_object()(v.cbegin(), v.cend()); + FT r = kernel_.squared_distance_d_object()(c, v[0]); + cache_.emplace_back(std::move(c), std::move(r)); + } + return cache_[k]; + } + + /// Return the circumradius, either from the old cache or computed, without writing to the cache. + template<class SimplicialComplexForAlpha> + auto radius(SimplicialComplexForAlpha& cplx, typename SimplicialComplexForAlpha::Simplex_handle s) { + auto k = cplx.key(s); + if(k!=cplx.null_key()) + return old_cache_[k].second; + // Using a transform_range is slower, currently. + thread_local std::vector<Point_d> v; + v.clear(); + for (auto vertex : cplx.simplex_vertex_range(s)) + v.push_back(get_point_(vertex)); + return kernel_.compute_squared_radius_d_object()(v.cbegin(), v.cend()); + } + public: /** \brief Inserts all Delaunay triangulation into the simplicial complex. - * It also computes the filtration values accordingly to the \ref createcomplexalgorithm + * It also computes the filtration values accordingly to the \ref createcomplexalgorithm if default_filtration_value + * is not set. * * \tparam SimplicialComplexForAlpha must meet `SimplicialComplexForAlpha` concept. * * @param[in] complex SimplicialComplexForAlpha to be created. * @param[in] max_alpha_square maximum for alpha square value. Default value is +\f$\infty\f$, and there is very - * little point using anything else since it does not save time. + * little point using anything else since it does not save time. Useless if `default_filtration_value` is set to + * `true`. * @param[in] exact Exact filtration values computation. Not exact if `Kernel` is not <a target="_blank" * href="https://doc.cgal.org/latest/Kernel_d/structCGAL_1_1Epeck__d.html">CGAL::Epeck_d</a>. - * + * @param[in] default_filtration_value Set this value to `true` if filtration values are not needed to be computed + * (will be set to `NaN`). + * Default value is `false` (which means compute the filtration values). + * * @return true if creation succeeds, false otherwise. * * @pre Delaunay triangulation must be already constructed with dimension strictly greater than 0. @@ -269,7 +319,8 @@ class Alpha_complex { typename Filtration_value = typename SimplicialComplexForAlpha::Filtration_value> bool create_complex(SimplicialComplexForAlpha& complex, Filtration_value max_alpha_square = std::numeric_limits<Filtration_value>::infinity(), - bool exact = false) { + bool exact = false, + bool default_filtration_value = false) { // From SimplicialComplexForAlpha type required to insert into a simplicial complex (with or without subfaces). typedef typename SimplicialComplexForAlpha::Vertex_handle Vertex_handle; typedef typename SimplicialComplexForAlpha::Simplex_handle Simplex_handle; @@ -296,19 +347,19 @@ class Alpha_complex { ++cit) { Vector_vertex vertexVector; #ifdef DEBUG_TRACES - std::cout << "Simplex_tree insertion "; + std::clog << "Simplex_tree insertion "; #endif // DEBUG_TRACES for (auto vit = cit->vertices_begin(); vit != cit->vertices_end(); ++vit) { if (*vit != nullptr) { #ifdef DEBUG_TRACES - std::cout << " " << (*vit)->data(); + std::clog << " " << (*vit)->data(); #endif // DEBUG_TRACES // Vector of vertex construction for simplex_tree structure vertexVector.push_back((*vit)->data()); } } #ifdef DEBUG_TRACES - std::cout << std::endl; + std::clog << std::endl; #endif // DEBUG_TRACES // Insert each simplex and its subfaces in the simplex tree - filtration is NaN complex.insert_simplex_and_subfaces(vertexVector, std::numeric_limits<double>::quiet_NaN()); @@ -316,62 +367,48 @@ class Alpha_complex { } // -------------------------------------------------------------------------------------------- - // -------------------------------------------------------------------------------------------- - // Will be re-used many times - Vector_of_CGAL_points pointVector; - // ### For i : d -> 0 - for (int decr_dim = triangulation_->maximal_dimension(); decr_dim >= 0; decr_dim--) { - // ### Foreach Sigma of dim i - for (Simplex_handle f_simplex : complex.skeleton_simplex_range(decr_dim)) { - int f_simplex_dim = complex.dimension(f_simplex); - if (decr_dim == f_simplex_dim) { - pointVector.clear(); -#ifdef DEBUG_TRACES - std::cout << "Sigma of dim " << decr_dim << " is"; -#endif // DEBUG_TRACES - for (auto vertex : complex.simplex_vertex_range(f_simplex)) { - pointVector.push_back(get_point(vertex)); -#ifdef DEBUG_TRACES - std::cout << " " << vertex; -#endif // DEBUG_TRACES - } -#ifdef DEBUG_TRACES - std::cout << std::endl; -#endif // DEBUG_TRACES - // ### If filt(Sigma) is NaN : filt(Sigma) = alpha(Sigma) - if (std::isnan(complex.filtration(f_simplex))) { - Filtration_value alpha_complex_filtration = 0.0; - // No need to compute squared_radius on a single point - alpha is 0.0 - if (f_simplex_dim > 0) { - // squared_radius function initialization - Squared_Radius squared_radius = kernel_.compute_squared_radius_d_object(); - - CGAL::NT_converter<typename Geom_traits::FT, Filtration_value> cv; - auto sqrad = squared_radius(pointVector.begin(), pointVector.end()); + if (!default_filtration_value) { + // -------------------------------------------------------------------------------------------- + // ### For i : d -> 0 + for (int decr_dim = triangulation_->maximal_dimension(); decr_dim >= 0; decr_dim--) { + // ### Foreach Sigma of dim i + for (Simplex_handle f_simplex : complex.skeleton_simplex_range(decr_dim)) { + int f_simplex_dim = complex.dimension(f_simplex); + if (decr_dim == f_simplex_dim) { + // ### If filt(Sigma) is NaN : filt(Sigma) = alpha(Sigma) + if (std::isnan(complex.filtration(f_simplex))) { + Filtration_value alpha_complex_filtration = 0.0; + // No need to compute squared_radius on a single point - alpha is 0.0 + if (f_simplex_dim > 0) { + auto const& sqrad = radius(complex, f_simplex); #if CGAL_VERSION_NR >= 1050000000 - if(exact) CGAL::exact(sqrad); + if(exact) CGAL::exact(sqrad); #endif - alpha_complex_filtration = cv(sqrad); - } - complex.assign_filtration(f_simplex, alpha_complex_filtration); + CGAL::NT_converter<FT, Filtration_value> cv; + alpha_complex_filtration = cv(sqrad); + } + complex.assign_filtration(f_simplex, alpha_complex_filtration); #ifdef DEBUG_TRACES - std::cout << "filt(Sigma) is NaN : filt(Sigma) =" << complex.filtration(f_simplex) << std::endl; + std::clog << "filt(Sigma) is NaN : filt(Sigma) =" << complex.filtration(f_simplex) << std::endl; #endif // DEBUG_TRACES + } + // No need to propagate further, unweighted points all have value 0 + if (decr_dim > 1) + propagate_alpha_filtration(complex, f_simplex); } - // No need to propagate further, unweighted points all have value 0 - if (decr_dim > 1) - propagate_alpha_filtration(complex, f_simplex); } + old_cache_ = std::move(cache_); + cache_.clear(); } + // -------------------------------------------------------------------------------------------- + + // -------------------------------------------------------------------------------------------- + // As Alpha value is an approximation, we have to make filtration non decreasing while increasing the dimension + complex.make_filtration_non_decreasing(); + // Remove all simplices that have a filtration value greater than max_alpha_square + complex.prune_above_filtration(max_alpha_square); + // -------------------------------------------------------------------------------------------- } - // -------------------------------------------------------------------------------------------- - - // -------------------------------------------------------------------------------------------- - // As Alpha value is an approximation, we have to make filtration non decreasing while increasing the dimension - complex.make_filtration_non_decreasing(); - // Remove all simplices that have a filtration value greater than max_alpha_square - complex.prune_above_filtration(max_alpha_square); - // -------------------------------------------------------------------------------------------- return true; } @@ -380,20 +417,18 @@ class Alpha_complex { void propagate_alpha_filtration(SimplicialComplexForAlpha& complex, Simplex_handle f_simplex) { // From SimplicialComplexForAlpha type required to assign filtration values. typedef typename SimplicialComplexForAlpha::Filtration_value Filtration_value; -#ifdef DEBUG_TRACES typedef typename SimplicialComplexForAlpha::Vertex_handle Vertex_handle; -#endif // DEBUG_TRACES // ### Foreach Tau face of Sigma for (auto f_boundary : complex.boundary_simplex_range(f_simplex)) { #ifdef DEBUG_TRACES - std::cout << " | --------------------------------------------------\n"; - std::cout << " | Tau "; + std::clog << " | --------------------------------------------------\n"; + std::clog << " | Tau "; for (auto vertex : complex.simplex_vertex_range(f_boundary)) { - std::cout << vertex << " "; + std::clog << vertex << " "; } - std::cout << "is a face of Sigma\n"; - std::cout << " | isnan(complex.filtration(Tau)=" << std::isnan(complex.filtration(f_boundary)) << std::endl; + std::clog << "is a face of Sigma\n"; + std::clog << " | isnan(complex.filtration(Tau)=" << std::isnan(complex.filtration(f_boundary)) << std::endl; #endif // DEBUG_TRACES // ### If filt(Tau) is not NaN if (!std::isnan(complex.filtration(f_boundary))) { @@ -402,37 +437,22 @@ class Alpha_complex { complex.filtration(f_simplex)); complex.assign_filtration(f_boundary, alpha_complex_filtration); #ifdef DEBUG_TRACES - std::cout << " | filt(Tau) = fmin(filt(Tau), filt(Sigma)) = " << complex.filtration(f_boundary) << std::endl; + std::clog << " | filt(Tau) = fmin(filt(Tau), filt(Sigma)) = " << complex.filtration(f_boundary) << std::endl; #endif // DEBUG_TRACES // ### Else } else { - // insert the Tau points in a vector for is_gabriel function - Vector_of_CGAL_points pointVector; -#ifdef DEBUG_TRACES - Vertex_handle vertexForGabriel = Vertex_handle(); -#endif // DEBUG_TRACES - for (auto vertex : complex.simplex_vertex_range(f_boundary)) { - pointVector.push_back(get_point(vertex)); - } - // Retrieve the Sigma point that is not part of Tau - parameter for is_gabriel function - Point_d point_for_gabriel; - for (auto vertex : complex.simplex_vertex_range(f_simplex)) { - point_for_gabriel = get_point(vertex); - if (std::find(pointVector.begin(), pointVector.end(), point_for_gabriel) == pointVector.end()) { -#ifdef DEBUG_TRACES - // vertex is not found in Tau - vertexForGabriel = vertex; -#endif // DEBUG_TRACES - // No need to continue loop - break; - } - } - // is_gabriel function initialization - Is_Gabriel is_gabriel = kernel_.side_of_bounded_sphere_d_object(); - bool is_gab = is_gabriel(pointVector.begin(), pointVector.end(), point_for_gabriel) - != CGAL::ON_BOUNDED_SIDE; + // Find which vertex of f_simplex is missing in f_boundary. We could actually write a variant of boundary_simplex_range that gives pairs (f_boundary, vertex). We rely on the fact that simplex_vertex_range is sorted. + auto longlist = complex.simplex_vertex_range(f_simplex); + auto shortlist = complex.simplex_vertex_range(f_boundary); + auto longiter = std::begin(longlist); + auto shortiter = std::begin(shortlist); + auto enditer = std::end(shortlist); + while(shortiter != enditer && *longiter == *shortiter) { ++longiter; ++shortiter; } + Vertex_handle extra = *longiter; + auto const& cache=get_cache(complex, f_boundary); + bool is_gab = kernel_.squared_distance_d_object()(cache.first, get_point_(extra)) >= cache.second; #ifdef DEBUG_TRACES - std::cout << " | Tau is_gabriel(Sigma)=" << is_gab << " - vertexForGabriel=" << vertexForGabriel << std::endl; + std::clog << " | Tau is_gabriel(Sigma)=" << is_gab << " - vertexForGabriel=" << extra << std::endl; #endif // DEBUG_TRACES // ### If Tau is not Gabriel of Sigma if (false == is_gab) { @@ -440,7 +460,7 @@ class Alpha_complex { Filtration_value alpha_complex_filtration = complex.filtration(f_simplex); complex.assign_filtration(f_boundary, alpha_complex_filtration); #ifdef DEBUG_TRACES - std::cout << " | filt(Tau) = filt(Sigma) = " << complex.filtration(f_boundary) << std::endl; + std::clog << " | filt(Tau) = filt(Sigma) = " << complex.filtration(f_boundary) << std::endl; #endif // DEBUG_TRACES } } diff --git a/src/Alpha_complex/include/gudhi/Alpha_complex_3d.h b/src/Alpha_complex/include/gudhi/Alpha_complex_3d.h index 7f96c94c..c19ebb79 100644 --- a/src/Alpha_complex/include/gudhi/Alpha_complex_3d.h +++ b/src/Alpha_complex/include/gudhi/Alpha_complex_3d.h @@ -61,10 +61,7 @@ namespace Gudhi { namespace alpha_complex { -#ifdef GUDHI_CAN_USE_CXX11_THREAD_LOCAL -thread_local -#endif // GUDHI_CAN_USE_CXX11_THREAD_LOCAL - double RELATIVE_PRECISION_OF_TO_DOUBLE = 0.00001; +thread_local double RELATIVE_PRECISION_OF_TO_DOUBLE = 0.00001; // Value_from_iterator returns the filtration value from an iterator on alpha shapes values // @@ -472,7 +469,7 @@ Weighted_alpha_complex_3d::Weighted_point_3 wp0(Weighted_alpha_complex_3d::Bare_ alpha_shape_3_ptr_->filtration_with_alpha_values(dispatcher); #ifdef DEBUG_TRACES - std::cout << "filtration_with_alpha_values returns : " << objects.size() << " objects" << std::endl; + std::clog << "filtration_with_alpha_values returns : " << objects.size() << " objects" << std::endl; #endif // DEBUG_TRACES using Alpha_value_iterator = typename std::vector<FT>::const_iterator; @@ -484,7 +481,7 @@ Weighted_alpha_complex_3d::Weighted_point_3 wp0(Weighted_alpha_complex_3d::Bare_ if (const Cell_handle* cell = CGAL::object_cast<Cell_handle>(&object_iterator)) { for (auto i = 0; i < 4; i++) { #ifdef DEBUG_TRACES - std::cout << "from cell[" << i << "] - Point coordinates (" << (*cell)->vertex(i)->point() << ")" + std::clog << "from cell[" << i << "] - Point coordinates (" << (*cell)->vertex(i)->point() << ")" << std::endl; #endif // DEBUG_TRACES vertex_list.push_back((*cell)->vertex(i)); @@ -496,7 +493,7 @@ Weighted_alpha_complex_3d::Weighted_point_3 wp0(Weighted_alpha_complex_3d::Bare_ for (auto i = 0; i < 4; i++) { if ((*facet).second != i) { #ifdef DEBUG_TRACES - std::cout << "from facet=[" << i << "] - Point coordinates (" << (*facet).first->vertex(i)->point() << ")" + std::clog << "from facet=[" << i << "] - Point coordinates (" << (*facet).first->vertex(i)->point() << ")" << std::endl; #endif // DEBUG_TRACES vertex_list.push_back((*facet).first->vertex(i)); @@ -508,7 +505,7 @@ Weighted_alpha_complex_3d::Weighted_point_3 wp0(Weighted_alpha_complex_3d::Bare_ } else if (const Edge* edge = CGAL::object_cast<Edge>(&object_iterator)) { for (auto i : {(*edge).second, (*edge).third}) { #ifdef DEBUG_TRACES - std::cout << "from edge[" << i << "] - Point coordinates (" << (*edge).first->vertex(i)->point() << ")" + std::clog << "from edge[" << i << "] - Point coordinates (" << (*edge).first->vertex(i)->point() << ")" << std::endl; #endif // DEBUG_TRACES vertex_list.push_back((*edge).first->vertex(i)); @@ -519,7 +516,7 @@ Weighted_alpha_complex_3d::Weighted_point_3 wp0(Weighted_alpha_complex_3d::Bare_ } else if (const Alpha_vertex_handle* vertex = CGAL::object_cast<Alpha_vertex_handle>(&object_iterator)) { #ifdef DEBUG_TRACES count_vertices++; - std::cout << "from vertex - Point coordinates (" << (*vertex)->point() << ")" << std::endl; + std::clog << "from vertex - Point coordinates (" << (*vertex)->point() << ")" << std::endl; #endif // DEBUG_TRACES vertex_list.push_back((*vertex)); } @@ -531,7 +528,7 @@ Weighted_alpha_complex_3d::Weighted_point_3 wp0(Weighted_alpha_complex_3d::Bare_ // alpha shape not found Complex_vertex_handle vertex = map_cgal_simplex_tree.size(); #ifdef DEBUG_TRACES - std::cout << "Point (" << the_alpha_shape_vertex->point() << ") not found - insert new vertex id " << vertex + std::clog << "Point (" << the_alpha_shape_vertex->point() << ") not found - insert new vertex id " << vertex << std::endl; #endif // DEBUG_TRACES the_simplex.push_back(vertex); @@ -540,7 +537,7 @@ Weighted_alpha_complex_3d::Weighted_point_3 wp0(Weighted_alpha_complex_3d::Bare_ // alpha shape found Complex_vertex_handle vertex = the_map_iterator->second; #ifdef DEBUG_TRACES - std::cout << "Point (" << the_alpha_shape_vertex->point() << ") found as vertex id " << vertex << std::endl; + std::clog << "Point (" << the_alpha_shape_vertex->point() << ") found as vertex id " << vertex << std::endl; #endif // DEBUG_TRACES the_simplex.push_back(vertex); } @@ -549,7 +546,7 @@ Weighted_alpha_complex_3d::Weighted_point_3 wp0(Weighted_alpha_complex_3d::Bare_ Filtration_value filtr = Value_from_iterator<Complexity>::perform(alpha_value_iterator); #ifdef DEBUG_TRACES - std::cout << "filtration = " << filtr << std::endl; + std::clog << "filtration = " << filtr << std::endl; #endif // DEBUG_TRACES complex.insert_simplex(the_simplex, static_cast<Filtration_value>(filtr)); GUDHI_CHECK(alpha_value_iterator != alpha_values.end(), "CGAL provided more simplices than values"); @@ -557,10 +554,10 @@ Weighted_alpha_complex_3d::Weighted_point_3 wp0(Weighted_alpha_complex_3d::Bare_ } #ifdef DEBUG_TRACES - std::cout << "vertices \t" << count_vertices << std::endl; - std::cout << "edges \t\t" << count_edges << std::endl; - std::cout << "facets \t\t" << count_facets << std::endl; - std::cout << "cells \t\t" << count_cells << std::endl; + std::clog << "vertices \t" << count_vertices << std::endl; + std::clog << "edges \t\t" << count_edges << std::endl; + std::clog << "facets \t\t" << count_facets << std::endl; + std::clog << "cells \t\t" << count_cells << std::endl; #endif // DEBUG_TRACES // -------------------------------------------------------------------------------------------- // As Alpha value is an approximation, we have to make filtration non decreasing while increasing the dimension diff --git a/src/Alpha_complex/test/Alpha_complex_3d_unit_test.cpp b/src/Alpha_complex/test/Alpha_complex_3d_unit_test.cpp index cd698a27..a4ecb6ad 100644 --- a/src/Alpha_complex/test/Alpha_complex_3d_unit_test.cpp +++ b/src/Alpha_complex/test/Alpha_complex_3d_unit_test.cpp @@ -54,7 +54,7 @@ BOOST_AUTO_TEST_CASE(Alpha_complex_3d_from_points) { // ----------------- // Fast version // ----------------- - std::cout << "Fast alpha complex 3d" << std::endl; + std::clog << "Fast alpha complex 3d" << std::endl; std::vector<Fast_alpha_complex_3d::Bare_point_3> points = get_points<Fast_alpha_complex_3d::Bare_point_3>(); Fast_alpha_complex_3d alpha_complex(points); @@ -79,7 +79,7 @@ BOOST_AUTO_TEST_CASE(Alpha_complex_3d_from_points) { // ----------------- // Exact version // ----------------- - std::cout << "Exact alpha complex 3d" << std::endl; + std::clog << "Exact alpha complex 3d" << std::endl; std::vector<Exact_alpha_complex_3d::Bare_point_3> exact_points = get_points<Exact_alpha_complex_3d::Bare_point_3>(); Exact_alpha_complex_3d exact_alpha_complex(exact_points); @@ -105,13 +105,13 @@ BOOST_AUTO_TEST_CASE(Alpha_complex_3d_from_points) { // --------------------- // Compare both versions // --------------------- - std::cout << "Exact Alpha complex 3d is of dimension " << exact_stree.dimension() << " - Fast is " + std::clog << "Exact Alpha complex 3d is of dimension " << exact_stree.dimension() << " - Fast is " << stree.dimension() << std::endl; BOOST_CHECK(exact_stree.dimension() == stree.dimension()); - std::cout << "Exact Alpha complex 3d num_simplices " << exact_stree.num_simplices() << " - Fast is " + std::clog << "Exact Alpha complex 3d num_simplices " << exact_stree.num_simplices() << " - Fast is " << stree.num_simplices() << std::endl; BOOST_CHECK(exact_stree.num_simplices() == stree.num_simplices()); - std::cout << "Exact Alpha complex 3d num_vertices " << exact_stree.num_vertices() << " - Fast is " + std::clog << "Exact Alpha complex 3d num_vertices " << exact_stree.num_vertices() << " - Fast is " << stree.num_vertices() << std::endl; BOOST_CHECK(exact_stree.num_vertices() == stree.num_vertices()); @@ -119,18 +119,18 @@ BOOST_AUTO_TEST_CASE(Alpha_complex_3d_from_points) { while (sh != stree.filtration_simplex_range().end()) { std::vector<int> simplex; std::vector<int> exact_simplex; - std::cout << "Fast ( "; + std::clog << "Fast ( "; for (auto vertex : stree.simplex_vertex_range(*sh)) { simplex.push_back(vertex); - std::cout << vertex << " "; + std::clog << vertex << " "; } - std::cout << ") -> [" << stree.filtration(*sh) << "] "; + std::clog << ") -> [" << stree.filtration(*sh) << "] "; // Find it in the exact structure auto sh_exact = exact_stree.find(simplex); BOOST_CHECK(sh_exact != exact_stree.null_simplex()); - std::cout << " versus [" << exact_stree.filtration(sh_exact) << "] " << std::endl; + std::clog << " versus [" << exact_stree.filtration(sh_exact) << "] " << std::endl; // Exact and non-exact version is not exactly the same due to float comparison GUDHI_TEST_FLOAT_EQUALITY_CHECK(exact_stree.filtration(sh_exact), stree.filtration(*sh)); @@ -139,7 +139,7 @@ BOOST_AUTO_TEST_CASE(Alpha_complex_3d_from_points) { // ----------------- // Safe version // ----------------- - std::cout << "Safe alpha complex 3d" << std::endl; + std::clog << "Safe alpha complex 3d" << std::endl; std::vector<Safe_alpha_complex_3d::Bare_point_3> safe_points = get_points<Safe_alpha_complex_3d::Bare_point_3>(); Safe_alpha_complex_3d safe_alpha_complex(safe_points); @@ -165,13 +165,13 @@ BOOST_AUTO_TEST_CASE(Alpha_complex_3d_from_points) { // --------------------- // Compare both versions // --------------------- - std::cout << "Safe Alpha complex 3d is of dimension " << safe_stree.dimension() << " - Fast is " + std::clog << "Safe Alpha complex 3d is of dimension " << safe_stree.dimension() << " - Fast is " << stree.dimension() << std::endl; BOOST_CHECK(safe_stree.dimension() == stree.dimension()); - std::cout << "Safe Alpha complex 3d num_simplices " << safe_stree.num_simplices() << " - Fast is " + std::clog << "Safe Alpha complex 3d num_simplices " << safe_stree.num_simplices() << " - Fast is " << stree.num_simplices() << std::endl; BOOST_CHECK(safe_stree.num_simplices() == stree.num_simplices()); - std::cout << "Safe Alpha complex 3d num_vertices " << safe_stree.num_vertices() << " - Fast is " + std::clog << "Safe Alpha complex 3d num_vertices " << safe_stree.num_vertices() << " - Fast is " << stree.num_vertices() << std::endl; BOOST_CHECK(safe_stree.num_vertices() == stree.num_vertices()); @@ -179,18 +179,18 @@ BOOST_AUTO_TEST_CASE(Alpha_complex_3d_from_points) { while (safe_sh != stree.filtration_simplex_range().end()) { std::vector<int> simplex; std::vector<int> exact_simplex; - std::cout << "Fast ( "; + std::clog << "Fast ( "; for (auto vertex : stree.simplex_vertex_range(*safe_sh)) { simplex.push_back(vertex); - std::cout << vertex << " "; + std::clog << vertex << " "; } - std::cout << ") -> [" << stree.filtration(*safe_sh) << "] "; + std::clog << ") -> [" << stree.filtration(*safe_sh) << "] "; // Find it in the exact structure auto sh_exact = safe_stree.find(simplex); BOOST_CHECK(sh_exact != safe_stree.null_simplex()); - std::cout << " versus [" << safe_stree.filtration(sh_exact) << "] " << std::endl; + std::clog << " versus [" << safe_stree.filtration(sh_exact) << "] " << std::endl; // Exact and non-exact version is not exactly the same due to float comparison GUDHI_TEST_FLOAT_EQUALITY_CHECK(safe_stree.filtration(sh_exact), stree.filtration(*safe_sh), 1e-15); diff --git a/src/Alpha_complex/test/Alpha_complex_unit_test.cpp b/src/Alpha_complex/test/Alpha_complex_unit_test.cpp index 27b671dd..4b37e4bd 100644 --- a/src/Alpha_complex/test/Alpha_complex_unit_test.cpp +++ b/src/Alpha_complex/test/Alpha_complex_unit_test.cpp @@ -48,7 +48,7 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(Alpha_complex_from_OFF_file, TestedKernel, list_of // ---------------------------------------------------------------------------- std::string off_file_name("alphacomplexdoc.off"); double max_alpha_square_value = 60.0; - std::cout << "========== OFF FILE NAME = " << off_file_name << " - alpha²=" << + std::clog << "========== OFF FILE NAME = " << off_file_name << " - alpha²=" << max_alpha_square_value << "==========" << std::endl; Gudhi::alpha_complex::Alpha_complex<TestedKernel> alpha_complex_from_file(off_file_name); @@ -56,29 +56,29 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(Alpha_complex_from_OFF_file, TestedKernel, list_of Gudhi::Simplex_tree<> simplex_tree_60; BOOST_CHECK(alpha_complex_from_file.create_complex(simplex_tree_60, max_alpha_square_value)); - std::cout << "simplex_tree_60.dimension()=" << simplex_tree_60.dimension() << std::endl; + std::clog << "simplex_tree_60.dimension()=" << simplex_tree_60.dimension() << std::endl; BOOST_CHECK(simplex_tree_60.dimension() == 2); - std::cout << "simplex_tree_60.num_vertices()=" << simplex_tree_60.num_vertices() << std::endl; + std::clog << "simplex_tree_60.num_vertices()=" << simplex_tree_60.num_vertices() << std::endl; BOOST_CHECK(simplex_tree_60.num_vertices() == 7); - std::cout << "simplex_tree_60.num_simplices()=" << simplex_tree_60.num_simplices() << std::endl; + std::clog << "simplex_tree_60.num_simplices()=" << simplex_tree_60.num_simplices() << std::endl; BOOST_CHECK(simplex_tree_60.num_simplices() == 25); max_alpha_square_value = 59.0; - std::cout << "========== OFF FILE NAME = " << off_file_name << " - alpha²=" << + std::clog << "========== OFF FILE NAME = " << off_file_name << " - alpha²=" << max_alpha_square_value << "==========" << std::endl; Gudhi::Simplex_tree<> simplex_tree_59; BOOST_CHECK(alpha_complex_from_file.create_complex(simplex_tree_59, max_alpha_square_value)); - std::cout << "simplex_tree_59.dimension()=" << simplex_tree_59.dimension() << std::endl; + std::clog << "simplex_tree_59.dimension()=" << simplex_tree_59.dimension() << std::endl; BOOST_CHECK(simplex_tree_59.dimension() == 2); - std::cout << "simplex_tree_59.num_vertices()=" << simplex_tree_59.num_vertices() << std::endl; + std::clog << "simplex_tree_59.num_vertices()=" << simplex_tree_59.num_vertices() << std::endl; BOOST_CHECK(simplex_tree_59.num_vertices() == 7); - std::cout << "simplex_tree_59.num_simplices()=" << simplex_tree_59.num_simplices() << std::endl; + std::clog << "simplex_tree_59.num_simplices()=" << simplex_tree_59.num_simplices() << std::endl; BOOST_CHECK(simplex_tree_59.num_simplices() == 23); } @@ -115,30 +115,30 @@ BOOST_AUTO_TEST_CASE(Alpha_complex_from_points) { // ---------------------------------------------------------------------------- Gudhi::alpha_complex::Alpha_complex<Kernel_4> alpha_complex_from_points(points); - std::cout << "========== Alpha_complex_from_points ==========" << std::endl; + std::clog << "========== Alpha_complex_from_points ==========" << std::endl; Gudhi::Simplex_tree<> simplex_tree; BOOST_CHECK(alpha_complex_from_points.create_complex(simplex_tree)); // Another way to check num_simplices - std::cout << "Iterator on alpha complex simplices in the filtration order, with [filtration value]:" << std::endl; + std::clog << "Iterator on alpha complex simplices in the filtration order, with [filtration value]:" << std::endl; int num_simplices = 0; for (auto f_simplex : simplex_tree.filtration_simplex_range()) { num_simplices++; - std::cout << " ( "; + std::clog << " ( "; for (auto vertex : simplex_tree.simplex_vertex_range(f_simplex)) { - std::cout << vertex << " "; + std::clog << vertex << " "; } - std::cout << ") -> " << "[" << simplex_tree.filtration(f_simplex) << "] "; - std::cout << std::endl; + std::clog << ") -> " << "[" << simplex_tree.filtration(f_simplex) << "] "; + std::clog << std::endl; } BOOST_CHECK(num_simplices == 15); - std::cout << "simplex_tree.num_simplices()=" << simplex_tree.num_simplices() << std::endl; + std::clog << "simplex_tree.num_simplices()=" << simplex_tree.num_simplices() << std::endl; BOOST_CHECK(simplex_tree.num_simplices() == 15); - std::cout << "simplex_tree.dimension()=" << simplex_tree.dimension() << std::endl; + std::clog << "simplex_tree.dimension()=" << simplex_tree.dimension() << std::endl; BOOST_CHECK(simplex_tree.dimension() == 3); - std::cout << "simplex_tree.num_vertices()=" << simplex_tree.num_vertices() << std::endl; + std::clog << "simplex_tree.num_vertices()=" << simplex_tree.num_vertices() << std::endl; BOOST_CHECK(simplex_tree.num_vertices() == points.size()); for (auto f_simplex : simplex_tree.filtration_simplex_range()) { @@ -162,22 +162,22 @@ BOOST_AUTO_TEST_CASE(Alpha_complex_from_points) { } Point_4 p0 = alpha_complex_from_points.get_point(0); - std::cout << "alpha_complex_from_points.get_point(0)=" << p0 << std::endl; + std::clog << "alpha_complex_from_points.get_point(0)=" << p0 << std::endl; BOOST_CHECK(4 == p0.dimension()); BOOST_CHECK(is_point_in_list(points, p0)); Point_4 p1 = alpha_complex_from_points.get_point(1); - std::cout << "alpha_complex_from_points.get_point(1)=" << p1 << std::endl; + std::clog << "alpha_complex_from_points.get_point(1)=" << p1 << std::endl; BOOST_CHECK(4 == p1.dimension()); BOOST_CHECK(is_point_in_list(points, p1)); Point_4 p2 = alpha_complex_from_points.get_point(2); - std::cout << "alpha_complex_from_points.get_point(2)=" << p2 << std::endl; + std::clog << "alpha_complex_from_points.get_point(2)=" << p2 << std::endl; BOOST_CHECK(4 == p2.dimension()); BOOST_CHECK(is_point_in_list(points, p2)); Point_4 p3 = alpha_complex_from_points.get_point(3); - std::cout << "alpha_complex_from_points.get_point(3)=" << p3 << std::endl; + std::clog << "alpha_complex_from_points.get_point(3)=" << p3 << std::endl; BOOST_CHECK(4 == p3.dimension()); BOOST_CHECK(is_point_in_list(points, p3)); @@ -188,30 +188,27 @@ BOOST_AUTO_TEST_CASE(Alpha_complex_from_points) { // Test after prune_above_filtration bool modified = simplex_tree.prune_above_filtration(0.6); - if (modified) { - simplex_tree.initialize_filtration(); - } BOOST_CHECK(modified); // Another way to check num_simplices - std::cout << "Iterator on alpha complex simplices in the filtration order, with [filtration value]:" << std::endl; + std::clog << "Iterator on alpha complex simplices in the filtration order, with [filtration value]:" << std::endl; num_simplices = 0; for (auto f_simplex : simplex_tree.filtration_simplex_range()) { num_simplices++; - std::cout << " ( "; + std::clog << " ( "; for (auto vertex : simplex_tree.simplex_vertex_range(f_simplex)) { - std::cout << vertex << " "; + std::clog << vertex << " "; } - std::cout << ") -> " << "[" << simplex_tree.filtration(f_simplex) << "] "; - std::cout << std::endl; + std::clog << ") -> " << "[" << simplex_tree.filtration(f_simplex) << "] "; + std::clog << std::endl; } BOOST_CHECK(num_simplices == 10); - std::cout << "simplex_tree.num_simplices()=" << simplex_tree.num_simplices() << std::endl; + std::clog << "simplex_tree.num_simplices()=" << simplex_tree.num_simplices() << std::endl; BOOST_CHECK(simplex_tree.num_simplices() == 10); - std::cout << "simplex_tree.dimension()=" << simplex_tree.dimension() << std::endl; + std::clog << "simplex_tree.dimension()=" << simplex_tree.dimension() << std::endl; BOOST_CHECK(simplex_tree.dimension() == 1); - std::cout << "simplex_tree.num_vertices()=" << simplex_tree.num_vertices() << std::endl; + std::clog << "simplex_tree.num_vertices()=" << simplex_tree.num_vertices() << std::endl; BOOST_CHECK(simplex_tree.num_vertices() == 4); for (auto f_simplex : simplex_tree.filtration_simplex_range()) { @@ -231,7 +228,7 @@ BOOST_AUTO_TEST_CASE(Alpha_complex_from_points) { } BOOST_AUTO_TEST_CASE_TEMPLATE(Alpha_complex_from_empty_points, TestedKernel, list_of_kernel_variants) { - std::cout << "========== Alpha_complex_from_empty_points ==========" << std::endl; + std::clog << "========== Alpha_complex_from_empty_points ==========" << std::endl; // ---------------------------------------------------------------------------- // Init of an empty list of points @@ -249,13 +246,13 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(Alpha_complex_from_empty_points, TestedKernel, lis Gudhi::Simplex_tree<> simplex_tree; BOOST_CHECK(!alpha_complex_from_points.create_complex(simplex_tree)); - std::cout << "simplex_tree.num_simplices()=" << simplex_tree.num_simplices() << std::endl; + std::clog << "simplex_tree.num_simplices()=" << simplex_tree.num_simplices() << std::endl; BOOST_CHECK(simplex_tree.num_simplices() == 0); - std::cout << "simplex_tree.dimension()=" << simplex_tree.dimension() << std::endl; + std::clog << "simplex_tree.dimension()=" << simplex_tree.dimension() << std::endl; BOOST_CHECK(simplex_tree.dimension() == -1); - std::cout << "simplex_tree.num_vertices()=" << simplex_tree.num_vertices() << std::endl; + std::clog << "simplex_tree.num_vertices()=" << simplex_tree.num_vertices() << std::endl; BOOST_CHECK(simplex_tree.num_vertices() == points.size()); } @@ -264,7 +261,7 @@ using Exact_kernel_2 = CGAL::Epeck_d< CGAL::Dimension_tag<2> >; using list_of_kernel_2_variants = boost::mpl::list<Inexact_kernel_2, Exact_kernel_2>; BOOST_AUTO_TEST_CASE_TEMPLATE(Alpha_complex_with_duplicated_points, TestedKernel, list_of_kernel_2_variants) { - std::cout << "========== Alpha_complex_with_duplicated_points ==========" << std::endl; + std::clog << "========== Alpha_complex_with_duplicated_points ==========" << std::endl; using Point = typename TestedKernel::Point_d; using Vector_of_points = std::vector<Point>; @@ -287,14 +284,14 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(Alpha_complex_with_duplicated_points, TestedKernel // ---------------------------------------------------------------------------- // Init of an alpha complex from the list of points // ---------------------------------------------------------------------------- - std::cout << "Init" << std::endl; + std::clog << "Init" << std::endl; Gudhi::alpha_complex::Alpha_complex<TestedKernel> alpha_complex_from_points(points); Gudhi::Simplex_tree<> simplex_tree; - std::cout << "create_complex" << std::endl; + std::clog << "create_complex" << std::endl; BOOST_CHECK(alpha_complex_from_points.create_complex(simplex_tree)); - std::cout << "simplex_tree.num_vertices()=" << simplex_tree.num_vertices() + std::clog << "simplex_tree.num_vertices()=" << simplex_tree.num_vertices() << std::endl; BOOST_CHECK(simplex_tree.num_vertices() < points.size()); } diff --git a/src/Alpha_complex/test/CMakeLists.txt b/src/Alpha_complex/test/CMakeLists.txt index 0476c6d4..fe4b23e4 100644 --- a/src/Alpha_complex/test/CMakeLists.txt +++ b/src/Alpha_complex/test/CMakeLists.txt @@ -8,11 +8,15 @@ if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0) add_executable ( Alpha_complex_test_unit Alpha_complex_unit_test.cpp ) target_link_libraries(Alpha_complex_test_unit ${CGAL_LIBRARY}) + add_executable ( Delaunay_complex_test_unit Delaunay_complex_unit_test.cpp ) + target_link_libraries(Delaunay_complex_test_unit ${CGAL_LIBRARY}) if (TBB_FOUND) target_link_libraries(Alpha_complex_test_unit ${TBB_LIBRARIES}) + target_link_libraries(Delaunay_complex_test_unit ${TBB_LIBRARIES}) endif() gudhi_add_boost_test(Alpha_complex_test_unit) + gudhi_add_boost_test(Delaunay_complex_test_unit) add_executable ( Alpha_complex_3d_test_unit Alpha_complex_3d_unit_test.cpp ) target_link_libraries(Alpha_complex_3d_test_unit ${CGAL_LIBRARY}) diff --git a/src/Alpha_complex/test/Delaunay_complex_unit_test.cpp b/src/Alpha_complex/test/Delaunay_complex_unit_test.cpp new file mode 100644 index 00000000..c1cc1fab --- /dev/null +++ b/src/Alpha_complex/test/Delaunay_complex_unit_test.cpp @@ -0,0 +1,68 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): Vincent Rouvreau + * + * Copyright (C) 2020 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#define BOOST_TEST_DYN_LINK +#define BOOST_TEST_MODULE "delaunay_complex" +#include <boost/test/unit_test.hpp> +#include <boost/mpl/list.hpp> + +#include <CGAL/Epick_d.h> +#include <CGAL/Epeck_d.h> + +#include <vector> +#include <limits> // NaN +#include <cmath> + +#include <gudhi/Alpha_complex.h> +// to construct a simplex_tree from Delaunay_triangulation +#include <gudhi/graph_simplicial_complex.h> +#include <gudhi/Simplex_tree.h> +#include <gudhi/Unitary_tests_utils.h> +#include <gudhi/random_point_generators.h> + +// Use dynamic_dimension_tag for the user to be able to set dimension +typedef CGAL::Epeck_d< CGAL::Dynamic_dimension_tag > Exact_kernel_d; +// Use static dimension_tag for the user not to be able to set dimension +typedef CGAL::Epeck_d< CGAL::Dimension_tag<5> > Exact_kernel_s; +// Use dynamic_dimension_tag for the user to be able to set dimension +typedef CGAL::Epick_d< CGAL::Dynamic_dimension_tag > Inexact_kernel_d; +// Use static dimension_tag for the user not to be able to set dimension +typedef CGAL::Epick_d< CGAL::Dimension_tag<5> > Inexact_kernel_s; +// The triangulation uses the default instantiation of the TriangulationDataStructure template parameter + +typedef boost::mpl::list<Exact_kernel_d, Exact_kernel_s, Inexact_kernel_d, Inexact_kernel_s> list_of_kernel_variants; + +using Simplex_tree = Gudhi::Simplex_tree<>; +using Simplex_handle = Simplex_tree::Simplex_handle; + +BOOST_AUTO_TEST_CASE_TEMPLATE(Alpha_complex_from_OFF_file, TestedKernel, list_of_kernel_variants) { + std::cout << "*****************************************************************************************************"; + using Point = typename TestedKernel::Point_d; + std::vector<Point> points; + // 50 points on a 4-sphere + points = Gudhi::generate_points_on_sphere_d<TestedKernel>(10, 5, 1.); + + Gudhi::alpha_complex::Alpha_complex<TestedKernel> alpha_complex(points); + + // Alpha complex + Simplex_tree stree_from_alpha_complex; + BOOST_CHECK(alpha_complex.create_complex(stree_from_alpha_complex)); + + // Delaunay complex + Simplex_tree stree_from_delaunay_complex; + BOOST_CHECK(alpha_complex.create_complex(stree_from_delaunay_complex, 0., false, true)); + + // Check all the simplices from alpha complex are in the Delaunay complex + for (auto f_simplex : stree_from_alpha_complex.complex_simplex_range()) { + Simplex_handle sh = stree_from_delaunay_complex.find(stree_from_alpha_complex.simplex_vertex_range(f_simplex)); + BOOST_CHECK(std::isnan(stree_from_delaunay_complex.filtration(sh))); + BOOST_CHECK(sh != stree_from_delaunay_complex.null_simplex()); + } +} diff --git a/src/Alpha_complex/test/Periodic_alpha_complex_3d_unit_test.cpp b/src/Alpha_complex/test/Periodic_alpha_complex_3d_unit_test.cpp index 731763fa..9eef920b 100644 --- a/src/Alpha_complex/test/Periodic_alpha_complex_3d_unit_test.cpp +++ b/src/Alpha_complex/test/Periodic_alpha_complex_3d_unit_test.cpp @@ -43,14 +43,14 @@ typedef boost::mpl::list<Fast_periodic_alpha_complex_3d, Safe_periodic_alpha_com periodic_variants_type_list; BOOST_AUTO_TEST_CASE_TEMPLATE(Alpha_complex_periodic_throw, Periodic_alpha_complex_3d, periodic_variants_type_list) { - std::cout << "Periodic alpha complex 3d exception throw" << std::endl; + std::clog << "Periodic alpha complex 3d exception throw" << std::endl; using Bare_point_3 = typename Periodic_alpha_complex_3d::Bare_point_3; std::vector<Bare_point_3> p_points; // Not important, this is not what we want to check p_points.push_back(Bare_point_3(0.0, 0.0, 0.0)); - std::cout << "Check exception throw in debug mode" << std::endl; + std::clog << "Check exception throw in debug mode" << std::endl; // Check it throws an exception when the cuboid is not iso BOOST_CHECK_THROW(Periodic_alpha_complex_3d periodic_alpha_complex(p_points, 0., 0., 0., 0.9, 1., 1.), std::invalid_argument); @@ -71,7 +71,7 @@ BOOST_AUTO_TEST_CASE(Alpha_complex_periodic) { // --------------------- // Fast periodic version // --------------------- - std::cout << "Fast periodic alpha complex 3d" << std::endl; + std::clog << "Fast periodic alpha complex 3d" << std::endl; using Creator = CGAL::Creator_uniform_3<double, Fast_periodic_alpha_complex_3d::Bare_point_3>; CGAL::Random random(7); @@ -106,7 +106,7 @@ BOOST_AUTO_TEST_CASE(Alpha_complex_periodic) { // ---------------------- // Exact periodic version // ---------------------- - std::cout << "Exact periodic alpha complex 3d" << std::endl; + std::clog << "Exact periodic alpha complex 3d" << std::endl; std::vector<Exact_periodic_alpha_complex_3d::Bare_point_3> e_p_points; @@ -122,13 +122,13 @@ BOOST_AUTO_TEST_CASE(Alpha_complex_periodic) { // --------------------- // Compare both versions // --------------------- - std::cout << "Exact periodic alpha complex 3d is of dimension " << exact_stree.dimension() << " - Non exact is " + std::clog << "Exact periodic alpha complex 3d is of dimension " << exact_stree.dimension() << " - Non exact is " << stree.dimension() << std::endl; BOOST_CHECK(exact_stree.dimension() == stree.dimension()); - std::cout << "Exact periodic alpha complex 3d num_simplices " << exact_stree.num_simplices() << " - Non exact is " + std::clog << "Exact periodic alpha complex 3d num_simplices " << exact_stree.num_simplices() << " - Non exact is " << stree.num_simplices() << std::endl; BOOST_CHECK(exact_stree.num_simplices() == stree.num_simplices()); - std::cout << "Exact periodic alpha complex 3d num_vertices " << exact_stree.num_vertices() << " - Non exact is " + std::clog << "Exact periodic alpha complex 3d num_vertices " << exact_stree.num_vertices() << " - Non exact is " << stree.num_vertices() << std::endl; BOOST_CHECK(exact_stree.num_vertices() == stree.num_vertices()); @@ -155,7 +155,7 @@ BOOST_AUTO_TEST_CASE(Alpha_complex_periodic) { // ---------------------- // Safe periodic version // ---------------------- - std::cout << "Safe periodic alpha complex 3d" << std::endl; + std::clog << "Safe periodic alpha complex 3d" << std::endl; std::vector<Safe_periodic_alpha_complex_3d::Bare_point_3> s_p_points; diff --git a/src/Alpha_complex/test/Weighted_alpha_complex_3d_unit_test.cpp b/src/Alpha_complex/test/Weighted_alpha_complex_3d_unit_test.cpp index 8035f6e8..6b31bea6 100644 --- a/src/Alpha_complex/test/Weighted_alpha_complex_3d_unit_test.cpp +++ b/src/Alpha_complex/test/Weighted_alpha_complex_3d_unit_test.cpp @@ -55,13 +55,13 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(Alpha_complex_weighted_throw, Weighted_alpha_compl // weights size is different from w_points size to make weighted Alpha_complex_3d throw in debug mode std::vector<double> weights = {0.01, 0.005, 0.006, 0.01, 0.009, 0.001}; - std::cout << "Check exception throw in debug mode" << std::endl; + std::clog << "Check exception throw in debug mode" << std::endl; BOOST_CHECK_THROW(Weighted_alpha_complex_3d wac(w_points, weights), std::invalid_argument); } #endif BOOST_AUTO_TEST_CASE_TEMPLATE(Alpha_complex_weighted, Weighted_alpha_complex_3d, weighted_variants_type_list) { - std::cout << "Weighted alpha complex 3d from points and weights" << std::endl; + std::clog << "Weighted alpha complex 3d from points and weights" << std::endl; using Bare_point_3 = typename Weighted_alpha_complex_3d::Bare_point_3; std::vector<Bare_point_3> w_points; w_points.push_back(Bare_point_3(0.0, 0.0, 0.0)); @@ -78,7 +78,7 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(Alpha_complex_weighted, Weighted_alpha_complex_3d, Gudhi::Simplex_tree<> stree; alpha_complex_p_a_w.create_complex(stree); - std::cout << "Weighted alpha complex 3d from weighted points" << std::endl; + std::clog << "Weighted alpha complex 3d from weighted points" << std::endl; using Weighted_point_3 = typename Weighted_alpha_complex_3d::Weighted_point_3; std::vector<Weighted_point_3> weighted_points; @@ -112,13 +112,13 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(Alpha_complex_weighted, Weighted_alpha_complex_3d, // --------------------- // Compare both versions // --------------------- - std::cout << "Weighted alpha complex 3d is of dimension " << stree_bis.dimension() << " - versus " + std::clog << "Weighted alpha complex 3d is of dimension " << stree_bis.dimension() << " - versus " << stree.dimension() << std::endl; BOOST_CHECK(stree_bis.dimension() == stree.dimension()); - std::cout << "Weighted alpha complex 3d num_simplices " << stree_bis.num_simplices() << " - versus " + std::clog << "Weighted alpha complex 3d num_simplices " << stree_bis.num_simplices() << " - versus " << stree.num_simplices() << std::endl; BOOST_CHECK(stree_bis.num_simplices() == stree.num_simplices()); - std::cout << "Weighted alpha complex 3d num_vertices " << stree_bis.num_vertices() << " - versus " + std::clog << "Weighted alpha complex 3d num_vertices " << stree_bis.num_vertices() << " - versus " << stree.num_vertices() << std::endl; BOOST_CHECK(stree_bis.num_vertices() == stree.num_vertices()); @@ -127,18 +127,18 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(Alpha_complex_weighted, Weighted_alpha_complex_3d, std::vector<int> simplex; std::vector<int> exact_simplex; #ifdef DEBUG_TRACES - std::cout << " ( "; + std::clog << " ( "; #endif for (auto vertex : stree.simplex_vertex_range(*sh)) { simplex.push_back(vertex); #ifdef DEBUG_TRACES - std::cout << vertex << " "; + std::clog << vertex << " "; #endif } #ifdef DEBUG_TRACES - std::cout << ") -> " + std::clog << ") -> " << "[" << stree.filtration(*sh) << "] "; - std::cout << std::endl; + std::clog << std::endl; #endif // Find it in the exact structure diff --git a/src/Alpha_complex/test/Weighted_periodic_alpha_complex_3d_unit_test.cpp b/src/Alpha_complex/test/Weighted_periodic_alpha_complex_3d_unit_test.cpp index b09e92d5..610b9f3d 100644 --- a/src/Alpha_complex/test/Weighted_periodic_alpha_complex_3d_unit_test.cpp +++ b/src/Alpha_complex/test/Weighted_periodic_alpha_complex_3d_unit_test.cpp @@ -45,7 +45,7 @@ typedef boost::mpl::list<Fast_weighted_periodic_alpha_complex_3d, Exact_weighted #ifdef GUDHI_DEBUG BOOST_AUTO_TEST_CASE_TEMPLATE(Alpha_complex_weighted_periodic_throw, Weighted_periodic_alpha_complex_3d, wp_variants_type_list) { - std::cout << "Weighted periodic alpha complex 3d exception throw" << std::endl; + std::clog << "Weighted periodic alpha complex 3d exception throw" << std::endl; using Creator = CGAL::Creator_uniform_3<double, typename Weighted_periodic_alpha_complex_3d::Bare_point_3>; CGAL::Random random(7); @@ -62,7 +62,7 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(Alpha_complex_weighted_periodic_throw, Weighted_pe p_weights.push_back(random.get_double(0., 0.01)); } - std::cout << "Cuboid is not iso exception" << std::endl; + std::clog << "Cuboid is not iso exception" << std::endl; // Check it throws an exception when the cuboid is not iso BOOST_CHECK_THROW( Weighted_periodic_alpha_complex_3d wp_alpha_complex(wp_points, p_weights, -1., -1., -1., 0.9, 1., 1.), @@ -83,7 +83,7 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(Alpha_complex_weighted_periodic_throw, Weighted_pe Weighted_periodic_alpha_complex_3d wp_alpha_complex(wp_points, p_weights, -1., -1., -1., 1., 1., 1.1), std::invalid_argument); - std::cout << "0 <= point.weight() < 1/64 * domain_size * domain_size exception" << std::endl; + std::clog << "0 <= point.weight() < 1/64 * domain_size * domain_size exception" << std::endl; // Weights must be in range ]0, 1/64 = 0.015625[ double temp = p_weights[25]; p_weights[25] = 1.0; @@ -97,7 +97,7 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(Alpha_complex_weighted_periodic_throw, Weighted_pe std::invalid_argument); p_weights[14] = temp; - std::cout << "wp_points and p_weights size exception" << std::endl; + std::clog << "wp_points and p_weights size exception" << std::endl; // Weights and points must have the same size // + 1 p_weights.push_back(1e-10); @@ -115,7 +115,7 @@ BOOST_AUTO_TEST_CASE(Alpha_complex_weighted_periodic) { // --------------------- // Fast weighted periodic version // --------------------- - std::cout << "Fast weighted periodic alpha complex 3d" << std::endl; + std::clog << "Fast weighted periodic alpha complex 3d" << std::endl; using Creator = CGAL::Creator_uniform_3<double, Fast_weighted_periodic_alpha_complex_3d::Bare_point_3>; CGAL::Random random(7); @@ -140,7 +140,7 @@ BOOST_AUTO_TEST_CASE(Alpha_complex_weighted_periodic) { // ---------------------- // Exact weighted periodic version // ---------------------- - std::cout << "Exact weighted periodic alpha complex 3d" << std::endl; + std::clog << "Exact weighted periodic alpha complex 3d" << std::endl; std::vector<Exact_weighted_periodic_alpha_complex_3d::Bare_point_3> e_p_points; @@ -156,13 +156,13 @@ BOOST_AUTO_TEST_CASE(Alpha_complex_weighted_periodic) { // --------------------- // Compare both versions // --------------------- - std::cout << "Exact weighted periodic alpha complex 3d is of dimension " << exact_stree.dimension() + std::clog << "Exact weighted periodic alpha complex 3d is of dimension " << exact_stree.dimension() << " - Non exact is " << stree.dimension() << std::endl; BOOST_CHECK(exact_stree.dimension() == stree.dimension()); - std::cout << "Exact weighted periodic alpha complex 3d num_simplices " << exact_stree.num_simplices() + std::clog << "Exact weighted periodic alpha complex 3d num_simplices " << exact_stree.num_simplices() << " - Non exact is " << stree.num_simplices() << std::endl; BOOST_CHECK(exact_stree.num_simplices() == stree.num_simplices()); - std::cout << "Exact weighted periodic alpha complex 3d num_vertices " << exact_stree.num_vertices() + std::clog << "Exact weighted periodic alpha complex 3d num_vertices " << exact_stree.num_vertices() << " - Non exact is " << stree.num_vertices() << std::endl; BOOST_CHECK(exact_stree.num_vertices() == stree.num_vertices()); @@ -189,7 +189,7 @@ BOOST_AUTO_TEST_CASE(Alpha_complex_weighted_periodic) { // ---------------------- // Safe weighted periodic version // ---------------------- - std::cout << "Safe weighted periodic alpha complex 3d" << std::endl; + std::clog << "Safe weighted periodic alpha complex 3d" << std::endl; std::vector<Safe_weighted_periodic_alpha_complex_3d::Bare_point_3> s_p_points; diff --git a/src/Alpha_complex/utilities/CMakeLists.txt b/src/Alpha_complex/utilities/CMakeLists.txt index 57b92942..2ffbdde0 100644 --- a/src/Alpha_complex/utilities/CMakeLists.txt +++ b/src/Alpha_complex/utilities/CMakeLists.txt @@ -2,7 +2,7 @@ project(Alpha_complex_utilities) if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0) add_executable (alpha_complex_persistence alpha_complex_persistence.cpp) - target_link_libraries(alpha_complex_persistence ${CGAL_LIBRARY} ${Boost_PROGRAM_OPTIONS_LIBRARY}) + target_link_libraries(alpha_complex_persistence ${CGAL_LIBRARY} Boost::program_options) if (TBB_FOUND) target_link_libraries(alpha_complex_persistence ${TBB_LIBRARIES}) @@ -16,14 +16,20 @@ if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0) if (DIFF_PATH) add_test(Alpha_complex_utilities_diff_exact_alpha_complex ${DIFF_PATH} "exact.pers" "safe.pers") + set_tests_properties(Alpha_complex_utilities_diff_exact_alpha_complex PROPERTIES DEPENDS + "Alpha_complex_utilities_exact_alpha_complex_persistence;Alpha_complex_utilities_safe_alpha_complex_persistence") + add_test(Alpha_complex_utilities_diff_fast_alpha_complex ${DIFF_PATH} "fast.pers" "safe.pers") + set_tests_properties(Alpha_complex_utilities_diff_fast_alpha_complex PROPERTIES DEPENDS + "Alpha_complex_utilities_fast_alpha_complex_persistence;Alpha_complex_utilities_safe_alpha_complex_persistence") + endif() install(TARGETS alpha_complex_persistence DESTINATION bin) add_executable(alpha_complex_3d_persistence alpha_complex_3d_persistence.cpp) - target_link_libraries(alpha_complex_3d_persistence ${CGAL_LIBRARY} ${Boost_PROGRAM_OPTIONS_LIBRARY}) + target_link_libraries(alpha_complex_3d_persistence ${CGAL_LIBRARY} Boost::program_options) if (TBB_FOUND) target_link_libraries(alpha_complex_3d_persistence ${TBB_LIBRARIES}) endif(TBB_FOUND) @@ -36,15 +42,19 @@ if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0) "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" "-p" "2" "-m" "0.45" "-o" "exact_3d.pers" "-e") - add_test(NAME Alpha_complex_utilities_safe_alpha_complex_3d COMMAND $<TARGET_FILE:alpha_complex_3d_persistence> + add_test(NAME Alpha_complex_utilities_fast_alpha_complex_3d COMMAND $<TARGET_FILE:alpha_complex_3d_persistence> "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" "-p" "2" "-m" "0.45" "-o" "fast_3d.pers" "-f") if (DIFF_PATH) add_test(Alpha_complex_utilities_diff_exact_alpha_complex_3d ${DIFF_PATH} "exact_3d.pers" "safe_3d.pers") + set_tests_properties(Alpha_complex_utilities_diff_exact_alpha_complex_3d PROPERTIES DEPENDS + "Alpha_complex_utilities_exact_alpha_complex_3d;Alpha_complex_utilities_alpha_complex_3d") add_test(Alpha_complex_utilities_diff_fast_alpha_complex_3d ${DIFF_PATH} "fast_3d.pers" "safe_3d.pers") + set_tests_properties(Alpha_complex_utilities_diff_fast_alpha_complex_3d PROPERTIES DEPENDS + "Alpha_complex_utilities_fast_alpha_complex_3d;Alpha_complex_utilities_alpha_complex_3d") endif() add_test(NAME Alpha_complex_utilities_periodic_alpha_complex_3d_persistence COMMAND $<TARGET_FILE:alpha_complex_3d_persistence> diff --git a/src/Alpha_complex/utilities/alpha_complex_3d_persistence.cpp b/src/Alpha_complex/utilities/alpha_complex_3d_persistence.cpp index 929fc2e8..91899040 100644 --- a/src/Alpha_complex/utilities/alpha_complex_3d_persistence.cpp +++ b/src/Alpha_complex/utilities/alpha_complex_3d_persistence.cpp @@ -222,10 +222,7 @@ int main(int argc, char **argv) { break; } - // Sort the simplices in the order of the filtration - simplex_tree.initialize_filtration(); - - std::cout << "Simplex_tree dim: " << simplex_tree.dimension() << std::endl; + std::clog << "Simplex_tree dim: " << simplex_tree.dimension() << std::endl; // Compute the persistence diagram of the complex Persistent_cohomology pcoh(simplex_tree, true); // initializes the coefficient field for homology @@ -237,7 +234,7 @@ int main(int argc, char **argv) { if (output_file_diag.empty()) { pcoh.output_diagram(); } else { - std::cout << "Result in file: " << output_file_diag << std::endl; + std::clog << "Result in file: " << output_file_diag << std::endl; std::ofstream out(output_file_diag); pcoh.output_diagram(out); out.close(); @@ -266,7 +263,7 @@ void program_options(int argc, char *argv[], std::string &off_file_points, bool "cuboid-file,c", po::value<std::string>(&cuboid_file), "Name of file describing the periodic domain. Format is:\n min_hx min_hy min_hz\n max_hx max_hy max_hz")( "output-file,o", po::value<std::string>(&output_file_diag)->default_value(std::string()), - "Name of file in which the persistence diagram is written. Default print in std::cout")( + "Name of file in which the persistence diagram is written. Default print in std::clog")( "max-alpha-square-value,r", po::value<Filtration_value>(&alpha_square_max_value) ->default_value(std::numeric_limits<Filtration_value>::infinity()), @@ -288,18 +285,18 @@ void program_options(int argc, char *argv[], std::string &off_file_points, bool po::notify(vm); if (vm.count("help") || !vm.count("input-file") || !vm.count("weight-file")) { - std::cout << std::endl; - std::cout << "Compute the persistent homology with coefficient field Z/pZ \n"; - std::cout << "of a 3D Alpha complex defined on a set of input points.\n"; - std::cout << "3D Alpha complex can be safe (by default) exact or fast, weighted and/or periodic\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.\n\n"; + std::clog << std::endl; + std::clog << "Compute the persistent homology with coefficient field Z/pZ \n"; + std::clog << "of a 3D Alpha complex defined on a set of input points.\n"; + std::clog << "3D Alpha complex can be safe (by default) exact or fast, weighted and/or periodic\n\n"; + std::clog << "The output diagram contains one bar per line, written with the convention: \n"; + std::clog << " p dim b d \n"; + std::clog << "where dim is the dimension of the homological feature,\n"; + std::clog << "b and d are respectively the birth and death of the feature and \n"; + std::clog << "p is the characteristic of the field Z/pZ used for homology coefficients.\n\n"; - std::cout << "Usage: " << argv[0] << " [options] input-file weight-file\n\n"; - std::cout << visible << std::endl; + std::clog << "Usage: " << argv[0] << " [options] input-file weight-file\n\n"; + std::clog << visible << std::endl; exit(-1); } } diff --git a/src/Alpha_complex/utilities/alpha_complex_persistence.cpp b/src/Alpha_complex/utilities/alpha_complex_persistence.cpp index 486347cc..7c898dfd 100644 --- a/src/Alpha_complex/utilities/alpha_complex_persistence.cpp +++ b/src/Alpha_complex/utilities/alpha_complex_persistence.cpp @@ -72,13 +72,10 @@ int main(int argc, char **argv) { // ---------------------------------------------------------------------------- // Display information about the alpha complex // ---------------------------------------------------------------------------- - std::cout << "Simplicial complex is of dimension " << simplex.dimension() << " - " << simplex.num_simplices() + std::clog << "Simplicial complex is of dimension " << simplex.dimension() << " - " << simplex.num_simplices() << " simplices - " << simplex.num_vertices() << " vertices." << std::endl; - // Sort the simplices in the order of the filtration - simplex.initialize_filtration(); - - std::cout << "Simplex_tree dim: " << simplex.dimension() << std::endl; + std::clog << "Simplex_tree dim: " << simplex.dimension() << std::endl; // Compute the persistence diagram of the complex Gudhi::persistent_cohomology::Persistent_cohomology<Simplex_tree, Gudhi::persistent_cohomology::Field_Zp> pcoh( simplex); @@ -91,7 +88,7 @@ int main(int argc, char **argv) { if (output_file_diag.empty()) { pcoh.output_diagram(); } else { - std::cout << "Result in file: " << output_file_diag << std::endl; + std::clog << "Result in file: " << output_file_diag << std::endl; std::ofstream out(output_file_diag); pcoh.output_diagram(out); out.close(); @@ -114,7 +111,7 @@ void program_options(int argc, char *argv[], std::string &off_file_points, bool "fast,f", po::bool_switch(&fast), "To activate fast version of Alpha complex (default is false, not available if exact is set)")( "output-file,o", po::value<std::string>(&output_file_diag)->default_value(std::string()), - "Name of file in which the persistence diagram is written. Default print in std::cout")( + "Name of file in which the persistence diagram is written. Default print in std::clog")( "max-alpha-square-value,r", po::value<Filtration_value>(&alpha_square_max_value) ->default_value(std::numeric_limits<Filtration_value>::infinity()), "Maximal alpha square value for the Alpha complex construction.")( @@ -135,17 +132,17 @@ void program_options(int argc, char *argv[], std::string &off_file_points, bool 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/pZ \n"; - std::cout << "of an Alpha 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::clog << std::endl; + std::clog << "Compute the persistent homology with coefficient field Z/pZ \n"; + std::clog << "of an Alpha complex defined on a set of input points.\n \n"; + std::clog << "The output diagram contains one bar per line, written with the convention: \n"; + std::clog << " p dim b d \n"; + std::clog << "where dim is the dimension of the homological feature,\n"; + std::clog << "b and d are respectively the birth and death of the feature and \n"; + std::clog << "p is the characteristic of the field Z/pZ used for homology coefficients." << std::endl << std::endl; + + std::clog << "Usage: " << argv[0] << " [options] input-file" << std::endl << std::endl; + std::clog << visible << std::endl; exit(-1); } } |