diff --git a/src/Alpha_complex/include/gudhi/Alpha_complex.h b/src/Alpha_complex/include/gudhi/Alpha_complex.h index ca9962046d..3074858db1 100644 --- a/src/Alpha_complex/include/gudhi/Alpha_complex.h +++ b/src/Alpha_complex/include/gudhi/Alpha_complex.h @@ -6,6 +6,7 @@ * * Modification(s): * - 2019/08 Vincent Rouvreau: Fix issue #10 for CGAL and Eigen3 + * - 2024/10 Vincent Rouvreau: Add square root filtration values interface * - YYYY/MM Author: Description of the modification */ @@ -65,14 +66,14 @@ template struct Is_Epeck_D> { static const bool val /** * \class Alpha_complex Alpha_complex.h gudhi/Alpha_complex.h * \brief Alpha complex data structure. - * + * * \ingroup alpha_complex - * + * * \details - * The data structure is constructing a CGAL Delaunay triangulation (for more information on CGAL Delaunay + * The data structure is constructing a CGAL Delaunay triangulation (for more information on CGAL Delaunay * triangulation, please refer to the corresponding chapter in page http://doc.cgal.org/latest/Triangulation/) from a * range of points or from an OFF file (cf. Points_off_reader). - * + * * Please refer to \ref alpha_complex for examples. * * The complex is a template class requiring an struct Is_Epeck_D> { static const bool val * href="https://doc.cgal.org/latest/Kernel_d/structCGAL_1_1Epeck__d.html">CGAL::Epeck_d * < * CGAL::Dynamic_dimension_tag > - * + * * \remark * - When Alpha_complex is constructed with an infinite value of alpha, the complex is a Delaunay complex. * - Using the default `CGAL::Epeck_d` makes the construction safe. If you pass exact=true to create_complex, the @@ -161,10 +162,10 @@ class Alpha_complex { public: /** \brief Alpha_complex constructor from an OFF file name. - * - * Uses the Points_off_reader to construct the Delaunay triangulation required to initialize + * + * Uses the Points_off_reader to construct the Delaunay triangulation required to initialize * the Alpha_complex. - * + * * Duplicate points are inserted once in the Alpha_complex. This is the reason why the vertices may be not contiguous. * * @param[in] off_file_name OFF file [path and] name. @@ -183,9 +184,9 @@ class Alpha_complex { * * The vertices may be not contiguous as some points may be discarded in the triangulation (duplicate points, * weighted hidden point, ...). - * + * * @param[in] points Range of points to triangulate. Points must be in Kernel::Point_d or Kernel::Weighted_point_d. - * + * * The type InputPointRange must be a range for which std::begin and std::end return input iterators on a * Kernel::Point_d or Kernel::Weighted_point_d. */ @@ -198,11 +199,11 @@ class Alpha_complex { * * The vertices may be not contiguous as some points may be discarded in the triangulation (duplicate points, * weighted hidden point, ...). - * + * * @param[in] points Range of points to triangulate. Points must be in Kernel::Point_d. - * + * * @param[in] weights Range of points weights. Weights must be in Kernel::FT. - * + * * The type InputPointRange must be a range for which std::begin and std::end return input iterators on a * Kernel::Point_d. */ @@ -354,8 +355,10 @@ class Alpha_complex { * It also computes the filtration values accordingly to the \ref createcomplexalgorithm if default_filtration_value * is not set. * + * \tparam square_root_filtrations If false (default value), it assigns to each simplex a filtration value equal to + * the squared cicumradius of the simplices, or to the radius when square_root_filtrations is true. * \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. Useless if `default_filtration_value` is set to @@ -367,13 +370,14 @@ class Alpha_complex { * 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. * @pre The simplicial complex must be empty (no vertices) - * + * * Initialization can be launched once. */ - template bool create_complex(SimplicialComplexForAlpha& complex, Filtration_value max_alpha_square = std::numeric_limits::infinity(), @@ -389,6 +393,9 @@ class Alpha_complex { using Simplex_handle = typename SimplicialComplexForAlpha::Simplex_handle; using Vector_vertex = std::vector; + // For users to be able to define their own sqrt function on their desired Filtration_value type + using std::sqrt; + if (triangulation_ == nullptr) { std::cerr << "Alpha_complex cannot create_complex from a NULL triangulation\n"; return false; // ----- >> @@ -438,7 +445,6 @@ class Alpha_complex { } } // -------------------------------------------------------------------------------------------- - if (!default_filtration_value) { CGAL::NT_converter cgal_converter; // -------------------------------------------------------------------------------------------- @@ -458,6 +464,9 @@ class Alpha_complex { if(exact) CGAL::exact(sqrad); #endif alpha_complex_filtration = cgal_converter(sqrad); + if constexpr (square_root_filtrations) { + alpha_complex_filtration = sqrt(alpha_complex_filtration); + } } complex.assign_filtration(f_simplex, alpha_complex_filtration); #ifdef DEBUG_TRACES @@ -473,14 +482,18 @@ class Alpha_complex { cache_.clear(); } // -------------------------------------------------------------------------------------------- - + // -------------------------------------------------------------------------------------------- if (!exact) // As Alpha value is an approximation, we have to make filtration non decreasing while increasing the dimension // Only in not exact version, cf. https://github.com/GUDHI/gudhi-devel/issues/57 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); + if constexpr (square_root_filtrations) { + complex.prune_above_filtration(sqrt(max_alpha_square)); + } else { + complex.prune_above_filtration(max_alpha_square); + } // -------------------------------------------------------------------------------------------- } return true; diff --git a/src/Alpha_complex/test/Delaunay_complex_unit_test.h b/src/Alpha_complex/test/Delaunay_complex_unit_test.h index e9a73f62da..df24441f6e 100644 --- a/src/Alpha_complex/test/Delaunay_complex_unit_test.h +++ b/src/Alpha_complex/test/Delaunay_complex_unit_test.h @@ -24,7 +24,7 @@ using Simplex_handle = Simplex_tree::Simplex_handle; template void compare_delaunay_complex_simplices() { - std::cout << "*****************************************************************************************************"; + std::clog << "*****************************************************************************************************\n"; using Point = typename CGAL_kernel::Point_d; std::vector points; // 50 points on a 4-sphere @@ -40,10 +40,24 @@ void compare_delaunay_complex_simplices() { 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 + std::clog << "Check all the simplices from alpha complex are in the Delaunay complex\n"; 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()); } + + std::clog << "Alpha complex with square root filtration values\n"; + Simplex_tree stree_from_alpha_sqrt; + BOOST_CHECK(alpha_complex.template create_complex(stree_from_alpha_sqrt)); + + std::clog << "Check simplices from alpha complex filtration values when square_root_filtrations is true\n"; + // 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_alpha_sqrt.find(stree_from_alpha_complex.simplex_vertex_range(f_simplex)); + BOOST_CHECK(sh != stree_from_alpha_sqrt.null_simplex()); + GUDHI_TEST_FLOAT_EQUALITY_CHECK(stree_from_alpha_sqrt.filtration(sh), + std::sqrt(stree_from_alpha_complex.filtration(f_simplex))); + } + } diff --git a/src/Alpha_complex/test/Weighted_alpha_complex_unit_test.cpp b/src/Alpha_complex/test/Weighted_alpha_complex_unit_test.cpp index 875704ee73..f40521761a 100644 --- a/src/Alpha_complex/test/Weighted_alpha_complex_unit_test.cpp +++ b/src/Alpha_complex/test/Weighted_alpha_complex_unit_test.cpp @@ -124,4 +124,47 @@ BOOST_AUTO_TEST_CASE(Weighted_alpha_complex_3d_comparison) { } ++dD_itr; } -} \ No newline at end of file +} + +BOOST_AUTO_TEST_CASE(Is_weighted_alpha_complex_nan) { + using Kernel = CGAL::Epeck_d< CGAL::Dimension_tag<3> >; + using Bare_point = Kernel::Point_d; + using Weighted_point = Kernel::Weighted_point_d; + using Vector_of_points = std::vector; + + Vector_of_points points; + 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.); + + Gudhi::alpha_complex::Alpha_complex alpha_complex_from_weighted_points(points); + + std::clog << "Weighted alpha complex\n"; + Gudhi::Simplex_tree<> stree; + if (alpha_complex_from_weighted_points.create_complex(stree)) { + for (auto f_simplex : stree.filtration_simplex_range()) { + std::clog << " ( "; + for (auto vertex : stree.simplex_vertex_range(f_simplex)) { + std::clog << vertex << " "; + } + std::clog << ") -> " << "[" << stree.filtration(f_simplex) << "]\n"; + + BOOST_CHECK(!std::isnan(stree.filtration(f_simplex))); + } + } + std::clog << "Weighted alpha complex with square_root_filtrations\n"; + Gudhi::Simplex_tree<> stree_sqrt; + if (alpha_complex_from_weighted_points.create_complex(stree_sqrt)) { + for (auto f_simplex : stree_sqrt.filtration_simplex_range()) { + std::clog << " ( "; + for (auto vertex : stree_sqrt.simplex_vertex_range(f_simplex)) { + std::clog << vertex << " "; + } + std::clog << ") -> " << "[" << stree_sqrt.filtration(f_simplex) << "]\n"; + + BOOST_CHECK(std::isnan(stree_sqrt.filtration(f_simplex))); + } + } +} diff --git a/src/Alpha_complex/utilities/CMakeLists.txt b/src/Alpha_complex/utilities/CMakeLists.txt index 237fcad2e1..519d68e772 100644 --- a/src/Alpha_complex/utilities/CMakeLists.txt +++ b/src/Alpha_complex/utilities/CMakeLists.txt @@ -7,6 +7,13 @@ if (TARGET CGAL::CGAL AND TARGET Eigen3::Eigen AND TARGET Boost::program_options "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" "-p" "2" "-m" "0.45" "-o" "fast.pers" "-f") add_test(NAME Alpha_complex_utilities_exact_alpha_complex_persistence COMMAND $ "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" "-p" "2" "-m" "0.45" "-o" "exact.pers" "-e") + # Same with sqrt filtration values - "-s" + add_test(NAME Alpha_complex_utilities_safe_sqrt_alpha_complex_persistence COMMAND $ + "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" "-s" "-p" "2" "-m" "0.45" "-o" "safe_sqrt.pers") + add_test(NAME Alpha_complex_utilities_fast_sqrt_alpha_complex_persistence COMMAND $ + "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" "-s" "-p" "2" "-m" "0.45" "-o" "fast_sqrt.pers" "-f") + add_test(NAME Alpha_complex_utilities_exact_sqrt_alpha_complex_persistence COMMAND $ + "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" "-s" "-p" "2" "-m" "0.45" "-o" "exact_sqrt.pers" "-e") if (DIFF_PATH) add_test(Alpha_complex_utilities_diff_exact_alpha_complex ${DIFF_PATH} "exact.pers" "safe.pers") @@ -17,6 +24,17 @@ if (TARGET CGAL::CGAL AND TARGET Eigen3::Eigen AND TARGET Boost::program_options "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") + + # Same with sqrt filtration values - "-s" + add_test(Alpha_complex_utilities_diff_exact_sqrt_alpha_complex ${DIFF_PATH} + "exact_sqrt.pers" "safe_sqrt.pers") + set_tests_properties(Alpha_complex_utilities_diff_exact_sqrt_alpha_complex PROPERTIES DEPENDS + "Alpha_complex_utilities_exact_sqrt_alpha_complex_persistence;Alpha_complex_utilities_safe_sqrt_alpha_complex_persistence") + + add_test(Alpha_complex_utilities_diff_fast_sqrt_alpha_complex ${DIFF_PATH} + "fast_sqrt.pers" "safe_sqrt.pers") + set_tests_properties(Alpha_complex_utilities_diff_fast_sqrt_alpha_complex PROPERTIES DEPENDS + "Alpha_complex_utilities_fast_sqrt_alpha_complex_persistence;Alpha_complex_utilities_safe_sqrt_alpha_complex_persistence") endif() install(TARGETS alpha_complex_persistence DESTINATION bin) @@ -62,6 +80,6 @@ if (TARGET CGAL::CGAL AND TARGET Boost::program_options) "-w" "${CMAKE_SOURCE_DIR}/data/points/grid_10_10_10_in_0_1.weights" "-c" "${CMAKE_SOURCE_DIR}/data/points/iso_cuboid_3_in_0_1.txt" "-p" "2" "-m" "0" "-e") - + install(TARGETS alpha_complex_3d_persistence DESTINATION bin) endif() diff --git a/src/Alpha_complex/utilities/alpha_complex_persistence.cpp b/src/Alpha_complex/utilities/alpha_complex_persistence.cpp index aad0ee9149..5cc4e997c5 100644 --- a/src/Alpha_complex/utilities/alpha_complex_persistence.cpp +++ b/src/Alpha_complex/utilities/alpha_complex_persistence.cpp @@ -5,6 +5,7 @@ * Copyright (C) 2016 Inria * * Modification(s): + * - 2024/10 Vincent Rouvreau: Add square root filtration values interface * - YYYY/MM Author: Description of the modification */ @@ -29,8 +30,9 @@ 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 &weight_file, std::string &output_file_diag, Filtration_value &alpha_square_max_value, - int &coeff_field_characteristic, Filtration_value &min_persistence); + bool& square_root_filtrations, 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) { @@ -61,10 +63,12 @@ std::vector read_weight_file(const std::string &weight_file) { 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) { + bool exact_version, Filtration_value alpha_square_max_value, + bool square_root_filtrations) { Simplex_tree stree; auto points = read_off(off_file_points); + bool complex_creation = false; if (weight_file != std::string()) { std::vector weights = read_weight_file(weight_file); if (points.size() != weights.size()) { @@ -75,18 +79,22 @@ Simplex_tree create_simplex_tree(const std::string &off_file_points, const std:: // 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); - } + if (square_root_filtrations) + complex_creation = alpha_complex_from_file.template create_complex(stree, alpha_square_max_value, exact_version); + else + complex_creation = alpha_complex_from_file.create_complex(stree, alpha_square_max_value, exact_version); } 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); - } + if (square_root_filtrations) + complex_creation = alpha_complex_from_file.template create_complex(stree, alpha_square_max_value, exact_version); + else + complex_creation = alpha_complex_from_file.create_complex(stree, alpha_square_max_value, exact_version); + } + if (!complex_creation) { + std::cerr << "Alpha complex simplicial complex creation failed." << std::endl; + exit(-1); } return stree; } @@ -97,12 +105,13 @@ int main(int argc, char **argv) { std::string output_file_diag; bool exact_version = false; bool fast_version = false; + bool square_root_filtrations = false; Filtration_value alpha_square_max_value; int coeff_field_characteristic; Filtration_value 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); + program_options(argc, argv, off_file_points, exact_version, fast_version, square_root_filtrations, 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; @@ -114,10 +123,10 @@ int main(int argc, char **argv) { // 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; - stree = create_simplex_tree(off_file_points, weight_file, exact_version, alpha_square_max_value); + stree = create_simplex_tree(off_file_points, weight_file, exact_version, alpha_square_max_value, square_root_filtrations); } else { using Kernel = CGAL::Epeck_d; - stree = create_simplex_tree(off_file_points, weight_file, exact_version, alpha_square_max_value); + stree = create_simplex_tree(off_file_points, weight_file, exact_version, alpha_square_max_value, square_root_filtrations); } // ---------------------------------------------------------------------------- // Display information about the alpha complex @@ -147,8 +156,9 @@ int main(int argc, char **argv) { } void program_options(int argc, char *argv[], std::string &off_file_points, bool &exact, bool &fast, - std::string &weight_file, std::string &output_file_diag, Filtration_value &alpha_square_max_value, - int &coeff_field_characteristic, Filtration_value &min_persistence) { + bool& square_root_filtrations, 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"); hidden.add_options()("input-file", po::value(&off_file_points), @@ -160,6 +170,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)")( + "square-root-filtrations,s", po::bool_switch(&square_root_filtrations), + "To activate square root filtration computations (default is false)")( "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()), @@ -191,7 +203,12 @@ void program_options(int argc, char *argv[], std::string &off_file_points, bool std::clog << " * fast: right combinatorics, values can be arbitrarily bad\n"; std::clog << " * safe (default): values can have a relative error at most 1e-5\n"; std::clog << " * exact: true values rounded to double.\n \n"; + std::clog << "Default Alpha complex filtrations computation are square of the circumradius of the simplex.\n"; + std::clog << "If you are interested in the circumradius of the simplex as filtration values, pass the "; + std::clog << "'-square-root-filtrations' (or '-s') option.\n"; std::clog << "Alpha complex can be, or not, weighted (requires a file containing weights values).\n\n"; + std::clog << "Weighted Alpha complex can have negative filtration values. If '-square-root-filtrations' is"; + std::clog << "set, filtration values will be Nan in this case.\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"; diff --git a/src/Alpha_complex/utilities/alphacomplex.md b/src/Alpha_complex/utilities/alphacomplex.md index eda9a469dd..5272b355e3 100644 --- a/src/Alpha_complex/utilities/alphacomplex.md +++ b/src/Alpha_complex/utilities/alphacomplex.md @@ -14,13 +14,19 @@ Leave the lines above as it is required by the web site generator 'Jekyll' This program computes the persistent homology with coefficient field *Z/pZ* of the dD alpha complex built from a dD point cloud. - + Different versions of Alpha complex computation are available: * fast: right combinatorics, values can be arbitrarily bad * safe (default): values can have a relative error at most 1e-5 * exact: true values rounded to double. - + +Default Alpha complex filtrations computation are square of the circumradius of the simplex. +If you are interested in the circumradius of the simplex as filtration values, pass the +'-square-root-filtrations' (or '-s') option. + Alpha complex can be, or not, weighted (requires a file containing weights values). +Weighted Alpha complex can have negative filtration values. If '-square-root-filtrations' is +set, filtration values will be Nan in this case. The output diagram contains one bar per line, written with the convention: @@ -82,7 +88,7 @@ Different versions of 3D Alpha complex computation are available: * fast: right combinatorics, values can be arbitrarily bad * safe (default): values can have a relative error at most 1e-5 * exact: true values rounded to double. - + 3D Alpha complex can be, or not, weighted (requires a file containing weights values) and/or periodic (requires a file describing the periodic domain). @@ -144,4 +150,4 @@ N.B.: and [Regular triangulation](https://doc.cgal.org/latest/Triangulation_3/index.html#Triangulation3secclassRegulartriangulation) documentation. * The periodic domain is detailed on CGAL [3D Periodic Triangulations User Manual]( -https://doc.cgal.org/latest/Periodic_3_triangulation_3/index.html) \ No newline at end of file +https://doc.cgal.org/latest/Periodic_3_triangulation_3/index.html) diff --git a/src/Cech_complex/include/gudhi/MEB_filtration.h b/src/Cech_complex/include/gudhi/MEB_filtration.h index 133c435fb1..cd30a7ead4 100644 --- a/src/Cech_complex/include/gudhi/MEB_filtration.h +++ b/src/Cech_complex/include/gudhi/MEB_filtration.h @@ -5,34 +5,44 @@ * Copyright (C) 2023 Inria * * Modification(s): + * - 2024/10 Vincent Rouvreau: Add square_root_filtrations argument to enable/disable squared radii computation * - YYYY/MM Author: Description of the modification */ #ifndef MEB_FILTRATION_H_ #define MEB_FILTRATION_H_ +#include + +#include +#include // for std::pair +#include // for std::sqrt + namespace Gudhi::cech_complex { /** * \ingroup cech_complex * * \brief - * Given a simplicial complex and an embedding of its vertices, this assigns to - * each simplex a filtration value equal to the squared radius of its minimal - * enclosing ball (MEB). + * Given a simplicial complex and an embedding of its vertices, this assigns to each simplex a filtration value equal + * to the squared (or not squared in function of `square_root_filtrations`) radius of its minimal enclosing ball (MEB). * - * Applied on a Čech complex, it recomputes the same values (squared). Applied on a Delaunay triangulation, it computes the Delaunay-Čech filtration. + * Applied on a Čech complex, it recomputes the same values (squared or not in function of `square_root_filtrations`). + * Applied on a Delaunay triangulation, it computes the Delaunay-Čech filtration. * + * \tparam square_root_filtrations If false (default value), it assigns to each simplex a filtration value equal to + * the squared radius of the MEB, or to the radius when square_root_filtrations is true. * \tparam Kernel CGAL kernel: either Epick_d or Epeck_d. * \tparam PointRange Random access range of `Kernel::Point_d`. * * @param[in] k The geometric kernel. * @param[in] complex The simplicial complex. * @param[in] points Embedding of the vertices of the complex. - * @param[in] exact If true and `Kernel` is CGAL::Epeck_d, the filtration values are computed exactly. Default is false. + * @param[in] exact If true and `Kernel` is + * CGAL::Epeck_d, the filtration values + * are computed exactly. Default is false. */ - -template +template void assign_MEB_filtration(Kernel&&k, SimplicialComplexForMEB& complex, PointRange const& points, bool exact = false) { using Point_d = typename Kernel::Point_d; using FT = typename Kernel::FT; @@ -42,6 +52,9 @@ void assign_MEB_filtration(Kernel&&k, SimplicialComplexForMEB& complex, PointRan using Simplex_handle = typename SimplicialComplexForMEB::Simplex_handle; using Filtration_value = typename SimplicialComplexForMEB::Filtration_value; + // For users to be able to define their own sqrt function on their desired Filtration_value type + using std::sqrt; + std::vector cache_; std::vector pts; CGAL::NT_converter cvt; @@ -68,7 +81,10 @@ void assign_MEB_filtration(Kernel&&k, SimplicialComplexForMEB& complex, PointRan FT r = k.squared_distance_d_object()(m, pu); if (exact) CGAL::exact(r); complex.assign_key(sh, cache_.size()); - complex.assign_filtration(sh, max(cvt(r), Filtration_value(0))); + Filtration_value filt{cvt(r)}; + if constexpr (square_root_filtrations) + filt = sqrt(filt); + complex.assign_filtration(sh, max(filt, Filtration_value(0))); cache_.emplace_back(std::move(m), std::move(r)); } else if (dim > ambient_dim) { // The sphere is always defined by at most d+1 points @@ -105,7 +121,11 @@ void assign_MEB_filtration(Kernel&&k, SimplicialComplexForMEB& complex, PointRan // int d2 = dim * dim; // Filtration_value max_sanity = maxf * d2 / (d2 - 1); // and use min(max_sanity, ...), which would limit how bad numerical errors can be. - maxf = max(maxf, cvt(r)); // maxf = cvt(r) except for rounding errors + Filtration_value filt{cvt(r)}; + if constexpr (square_root_filtrations) + filt = sqrt(filt); + // maxf = filt except for rounding errors + maxf = max(maxf, filt); complex.assign_key(sh, cache_.size()); // We could check if the simplex is maximal and avoiding adding it to the cache in that case. cache_.emplace_back(std::move(c), std::move(r)); diff --git a/src/Cech_complex/test/test_cech_complex.cpp b/src/Cech_complex/test/test_cech_complex.cpp index f5b6cbc804..8d0901834e 100644 --- a/src/Cech_complex/test/test_cech_complex.cpp +++ b/src/Cech_complex/test/test_cech_complex.cpp @@ -174,10 +174,34 @@ BOOST_AUTO_TEST_CASE(Cech_complex_for_documentation) { cech5.create_complex(st3, 5); BOOST_CHECK(st3.dimension() == 5); auto st3_save = st3; + std::clog << "Simplex tree from Cech complex\n"; + for (auto f_simplex : st3.filtration_simplex_range()) { + std::clog << " ( "; + for (auto vertex : st3.simplex_vertex_range(f_simplex)) std::clog << vertex << " "; + std::clog << ") -> " << "[" << st3.filtration(f_simplex) << "]\n"; + } + st3.reset_filtration(-1); // unnecessary, but ensures we don't cheat Gudhi::cech_complex::assign_MEB_filtration(Kernel(), st3, points); for (auto sh : st3.complex_simplex_range()) st3.assign_filtration(sh, std::sqrt(st3.filtration(sh))); + std::clog << "Simplex tree from assign_MEB_filtration - after std::sqrt on filtration values\n"; + for (auto f_simplex : st3.filtration_simplex_range()) { + std::clog << " ( "; + for (auto vertex : st3.simplex_vertex_range(f_simplex)) std::clog << vertex << " "; + std::clog << ") -> " << "[" << st3.filtration(f_simplex) << "]\n"; + } + BOOST_CHECK(st3 == st3_save); // Should only be an approximate test + + // Same test but with square_root_filtrations=true + st3.reset_filtration(-1); // unnecessary, but ensures we don't cheat + Gudhi::cech_complex::assign_MEB_filtration(Kernel(), st3, points); + std::clog << "Simplex tree from assign_MEB_filtration with square_root_filtrations=true\n"; + for (auto f_simplex : st3.filtration_simplex_range()) { + std::clog << " ( "; + for (auto vertex : st3.simplex_vertex_range(f_simplex)) std::clog << vertex << " "; + std::clog << ") -> " << "[" << st3.filtration(f_simplex) << "]\n"; + } BOOST_CHECK(st3 == st3_save); // Should only be an approximate test } diff --git a/src/Cech_complex/utilities/CMakeLists.txt b/src/Cech_complex/utilities/CMakeLists.txt index 44a1d49de8..1c8850545f 100644 --- a/src/Cech_complex/utilities/CMakeLists.txt +++ b/src/Cech_complex/utilities/CMakeLists.txt @@ -29,6 +29,13 @@ if (TARGET CGAL::CGAL AND TARGET Eigen3::Eigen AND TARGET Boost::program_options "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" "-r" "0.25" "-m" "0.5" "-p" "3" "-o" "delaunay_fast.pers" "-f") add_test(NAME Delaunay_Cech_complex_utility_from_rips_on_tore_3D_exact COMMAND $ "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" "-r" "0.25" "-m" "0.5" "-p" "3" "-o" "delaunay_exact.pers" "-e") + # Same with sqrt filtration values - "-s" + add_test(NAME Delaunay_Cech_complex_utility_from_rips_on_tore_3D_sqrt_safe COMMAND $ + "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" "-s" "-r" "0.25" "-m" "0.5" "-p" "3" "-o" "delaunay_sqrt_safe.pers") + add_test(NAME Delaunay_Cech_complex_utility_from_rips_on_tore_3D_sqrt_fast COMMAND $ + "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" "-s" "-r" "0.25" "-m" "0.5" "-p" "3" "-o" "delaunay_sqrt_fast.pers" "-f") + add_test(NAME Delaunay_Cech_complex_utility_from_rips_on_tore_3D_sqrt_exact COMMAND $ + "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" "-s" "-r" "0.25" "-m" "0.5" "-p" "3" "-o" "delaunay_sqrt_exact.pers" "-e") if (DIFF_PATH) add_test(Delaunay_Cech_complex_utilities_diff_exact ${DIFF_PATH} @@ -40,6 +47,17 @@ if (TARGET CGAL::CGAL AND TARGET Eigen3::Eigen AND TARGET Boost::program_options "delaunay_fast.pers" "delaunay_safe.pers") set_tests_properties(Delaunay_Cech_complex_utilities_diff_fast PROPERTIES DEPENDS "Delaunay_Cech_complex_utility_from_rips_on_tore_3D_safe;Delaunay_Cech_complex_utility_from_rips_on_tore_3D_fast") + + # Same with sqrt filtration values - "-s" + add_test(Delaunay_Cech_complex_utilities_diff_sqrt_exact ${DIFF_PATH} + "delaunay_sqrt_exact.pers" "delaunay_sqrt_safe.pers") + set_tests_properties(Delaunay_Cech_complex_utilities_diff_sqrt_exact PROPERTIES DEPENDS + "Delaunay_Cech_complex_utility_from_rips_on_tore_3D_sqrt_safe;Delaunay_Cech_complex_utility_from_rips_on_tore_3D_sqrt_exact") + + add_test(Delaunay_Cech_complex_utilities_diff_sqrt_fast ${DIFF_PATH} + "delaunay_sqrt_fast.pers" "delaunay_sqrt_safe.pers") + set_tests_properties(Delaunay_Cech_complex_utilities_diff_sqrt_fast PROPERTIES DEPENDS + "Delaunay_Cech_complex_utility_from_rips_on_tore_3D_sqrt_safe;Delaunay_Cech_complex_utility_from_rips_on_tore_3D_sqrt_fast") endif() install(TARGETS delaunay_cech_persistence DESTINATION bin) diff --git a/src/Cech_complex/utilities/cechcomplex.md b/src/Cech_complex/utilities/cechcomplex.md index 4201e0ed69..290650d9df 100644 --- a/src/Cech_complex/utilities/cechcomplex.md +++ b/src/Cech_complex/utilities/cechcomplex.md @@ -62,7 +62,7 @@ Beware: this program may use a lot of RAM and take a lot of time if `max-radius` ## dealaunay_cech_persistence ## -This program Computes the persistent homology with coefficient field *Z/pZ* +This program Computes the persistent homology with coefficient field *Z/pZ* of a Delaunay-Čech complex defined on a set of input points. Different versions of Delaunay-Čech complex computation are available: @@ -70,7 +70,13 @@ Different versions of Delaunay-Čech complex computation are available: * safe (default): values can have a relative error at most 1e-5 * exact: true values rounded to double. -The output diagram contains one bar per line, written with the + +Default Delaunay-Čech complex filtrations computation are squared radius of the MEB. +If you are interested in radius of the MEB as filtration values, pass the '-square-root-filtrations' +(or '-s') option. + + + The output diagram contains one bar per line, written with the convention: `p dim birth death` diff --git a/src/Cech_complex/utilities/delaunay_cech_persistence.cpp b/src/Cech_complex/utilities/delaunay_cech_persistence.cpp index d7417f12af..e967a2b5bb 100644 --- a/src/Cech_complex/utilities/delaunay_cech_persistence.cpp +++ b/src/Cech_complex/utilities/delaunay_cech_persistence.cpp @@ -5,6 +5,7 @@ * Copyright (C) 2024 Inria * * Modification(s): + * - 2024/10 Vincent Rouvreau: Add square_root_filtrations argument to enable/disable squared radii computation * - YYYY/MM Author: Description of the modification */ @@ -31,15 +32,18 @@ using Field_Zp = Gudhi::persistent_cohomology::Field_Zp; using Persistent_cohomology = Gudhi::persistent_cohomology::Persistent_cohomology; void program_options(int argc, char* argv[], std::string& off_file_points, bool& exact, bool& fast, - std::string& pers_diag_file, Filtration_value& max_radius, int& p, + bool& square_root_filtrations, std::string& pers_diag_file, Filtration_value& max_radius, int& p, Filtration_value& min_persistence); template -Simplex_tree create_simplex_tree(const std::string &off_file_points, bool exact_version, Filtration_value max_radius) { +Simplex_tree create_simplex_tree(const std::string &off_file_points, bool exact_version, Filtration_value max_radius, + bool square_root_filtrations) { using Point = typename Kernel::Point_d; using Points_off_reader = Gudhi::Points_off_reader; using Delaunay_complex = Gudhi::alpha_complex::Alpha_complex; + using std::sqrt; + Simplex_tree stree; Points_off_reader off_reader(off_file_points); @@ -48,7 +52,12 @@ Simplex_tree create_simplex_tree(const std::string &off_file_points, bool exact_ delaunay_complex_from_file.create_complex(stree, std::numeric_limits< Filtration_value >::infinity(), // exact can be false (or true), as default_filtration_value is set to true false, true); - Gudhi::cech_complex::assign_MEB_filtration(Kernel(), stree, point_cloud, exact_version); + if (square_root_filtrations) { + Gudhi::cech_complex::assign_MEB_filtration(Kernel(), stree, point_cloud, exact_version); + max_radius = sqrt(max_radius); + } else { + Gudhi::cech_complex::assign_MEB_filtration(Kernel(), stree, point_cloud, exact_version); + } stree.prune_above_filtration(max_radius); return stree; } @@ -58,12 +67,13 @@ int main(int argc, char* argv[]) { std::string pers_diag_file; bool exact_version = false; bool fast_version = false; + bool square_root_filtrations = false; Filtration_value max_radius; int p; Filtration_value min_persistence; - program_options(argc, argv, off_file_points, exact_version, fast_version, pers_diag_file, max_radius, p, - min_persistence); + program_options(argc, argv, off_file_points, exact_version, fast_version, square_root_filtrations, pers_diag_file, + max_radius, p, min_persistence); if ((exact_version) && (fast_version)) { std::cerr << "You cannot set the exact and the fast version." << std::endl; @@ -75,11 +85,11 @@ int main(int argc, char* argv[]) { // 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; - stree = create_simplex_tree(off_file_points, exact_version, max_radius); + stree = create_simplex_tree(off_file_points, exact_version, max_radius, square_root_filtrations); } else { std::clog << "exact_version = " << exact_version << "\n"; using Kernel = CGAL::Epeck_d; - stree = create_simplex_tree(off_file_points, exact_version, max_radius); + stree = create_simplex_tree(off_file_points, exact_version, max_radius, square_root_filtrations); } std::clog << "The complex contains " << stree.num_simplices() << " simplices \n"; @@ -108,7 +118,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& pers_diag_file, Filtration_value& max_radius, int& p, + bool& square_root_filtrations, std::string& pers_diag_file, Filtration_value& max_radius, int& p, Filtration_value& min_persistence) { namespace po = boost::program_options; po::options_description hidden("Hidden options"); @@ -121,6 +131,8 @@ void program_options(int argc, char* argv[], std::string& off_file_points, bool& "To activate exact version of Delaunay-Cech complex (default is false, not available if fast is set)")( "fast,f", po::bool_switch(&fast), "To activate fast version of Delaunay-Cech complex (default is false, not available if exact is set)")( + "square-root-filtrations,s", po::bool_switch(&square_root_filtrations), + "To activate square root filtration computations (default is false)")( "output-file,o", po::value(&pers_diag_file)->default_value(std::string()), "Name of file in which the persistence diagram is written. Default print in standard output")( "max-radius,r", @@ -150,6 +162,9 @@ void program_options(int argc, char* argv[], std::string& off_file_points, bool& std::clog << " * fast: right combinatorics, values can be arbitrarily bad\n"; std::clog << " * safe (default): values can have a relative error at most 1e-5\n"; std::clog << " * exact: true values rounded to double.\n \n"; + std::clog << "Default Delaunay-Cech complex filtrations computation are squared radius of the MEB.\n"; + std::clog << "If you are interested in radius of the MEB as filtration values, pass the '-square-root-filtrations'"; + std::clog << "(or '-s') option.\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"; diff --git a/src/python/gudhi/delaunay_complex.pyx b/src/python/gudhi/delaunay_complex.pyx index 49144c89fa..8938b57e6b 100644 --- a/src/python/gudhi/delaunay_complex.pyx +++ b/src/python/gudhi/delaunay_complex.pyx @@ -8,6 +8,7 @@ # # Modification(s): # - 2024/03 Vincent Rouvreau: Renamed AlphaComplex as DelaunayComplex. AlphaComplex inherits from it. +# - 2024/10 Vincent Rouvreau: Add square root filtration values interface # - YYYY/MM Author: Description of the modification from __future__ import print_function @@ -18,6 +19,7 @@ from libcpp.string cimport string from libcpp cimport bool from libc.stdint cimport intptr_t import warnings +from typing import Literal, Optional, Iterable from gudhi.simplex_tree cimport * from gudhi.simplex_tree import SimplexTree @@ -40,7 +42,7 @@ cdef extern from "Delaunay_complex_interface.h" namespace "Gudhi": cdef cppclass Delaunay_complex_interface "Gudhi::delaunay_complex::Delaunay_complex_interface": Delaunay_complex_interface(vector[vector[double]] points, vector[double] weights, bool fast_version, bool exact_version) nogil except + vector[double] get_point(int vertex) nogil except + - void create_simplex_tree(Simplex_tree_python_interface* simplex_tree, double max_alpha_square, Delaunay_filtration filtration) nogil except + + void create_simplex_tree(Simplex_tree_python_interface* simplex_tree, double max_alpha_square, Delaunay_filtration filtration, bool square_root_filtrations) nogil except + @staticmethod void set_float_relative_precision(double precision) nogil @staticmethod @@ -60,23 +62,21 @@ cdef class DelaunayComplex: cdef Delaunay_complex_interface * this_ptr # Fake constructor that does nothing but documenting the constructor - def __init__(self, points=[], weights=None, precision='safe'): + def __init__(self, points: Iterable[Iterable[float]] = [], weights: Optional[Iterable[float]] = None, + precision: Literal['fast', 'safe', 'exact'] = 'safe'): """DelaunayComplex constructor. - :param points: A list of points in d-Dimension. - :type points: Iterable[Iterable[float]] - - :param weights: A list of weights. If set, the number of weights must correspond to the number of points. - :type weights: Iterable[float] - - :param precision: Delaunay complex precision can be 'fast', 'safe' or 'exact'. Default is 'safe'. - :type precision: string + Args: + points (Iterable[Iterable[float]]): A list of points in d-Dimension. + weights (Optional[Iterable[float]]): A list of weights. If set, the number of weights must correspond to the number of points. + precision (str): Delaunay complex precision can be 'fast', 'safe' or 'exact'. Default is 'safe'. :raises ValueError: In case of inconsistency between the number of points and weights. """ # The real cython constructor - def __cinit__(self, points = [], weights=None, precision = 'safe'): + def __cinit__(self, points: Iterable[Iterable[float]] = [], weights: Optional[Iterable[float]] = None, + precision: Literal['fast', 'safe', 'exact'] = 'safe'): assert precision in ['fast', 'safe', 'exact'], "Delaunay complex precision can only be 'fast', 'safe' or 'exact'" cdef bool fast = precision == 'fast' cdef bool exact = precision == 'exact' @@ -103,25 +103,28 @@ cdef class DelaunayComplex: """ return self.this_ptr != NULL - def create_simplex_tree(self, max_alpha_square = float('inf'), filtration = None): + def create_simplex_tree(self, max_alpha_square: float = float('inf'), + filtration: Optional[Literal['alpha', 'cech']] = None, + square_root_filtrations: bool = False): """ - :param max_alpha_square: The maximum alpha square threshold the simplices shall not exceed. Default is set to - infinity, and there is very little point using anything else since it does not save time. - :type max_alpha_square: float - :param filtration: Set this value to `None` (default value) if filtration values are not needed to be computed - (will be set to `NaN`). Set it to `alpha` to compute the filtration values with the Alpha complex, or to - `cech` to compute the Delaunay Cech complex. - :type filtration: string or None - :returns: A simplex tree created from the Delaunay Triangulation. The vertex `k` corresponds to the k-th input - point. The vertices may not be numbered contiguously as some points may be discarded in the triangulation - (duplicate points, weighted hidden point, ...). - :rtype: SimplexTree + Args: + max_alpha_square: The maximum alpha square threshold the simplices shall not exceed. Default is set to + infinity, and there is very little point using anything else since it does not save time. + filtration: Set this value to `None` (default value) if filtration values are not needed to be computed + (will be set to `NaN`). Set it to `alpha` to compute the filtration values with the Alpha complex, or + to `cech` to compute the Delaunay Cech complex. + square_root_filtrations: Square root filtration values when True. Default is False. + Returns: + SimplexTree: A simplex tree created from the Delaunay Triangulation. The vertex `k` corresponds to the k-th input + point. The vertices may not be numbered contiguously as some points may be discarded in the + triangulation (duplicate points, weighted hidden point, ...). """ if not filtration in [None, 'alpha', 'cech']: raise ValueError(f"\'{filtration}\' is not a valid filtration value. Must be None, \'alpha\' or \'cech\'") stree = SimplexTree() cdef double mas = max_alpha_square cdef intptr_t stree_int_ptr=stree.thisptr + cdef bool srf = square_root_filtrations cdef Delaunay_filtration filt = NONE if filtration == 'cech': @@ -131,31 +134,31 @@ cdef class DelaunayComplex: with nogil: self.this_ptr.create_simplex_tree(stree_int_ptr, - mas, filt) + mas, filt, srf) return stree @staticmethod - def set_float_relative_precision(precision): + def set_float_relative_precision(precision: float): """ - :param precision: When the DelaunayComplex is constructed with :code:`precision = 'safe'` (the default), - one can set the float relative precision of filtration values computed in - :func:`~gudhi.DelaunayComplex.create_simplex_tree`. - Default is :code:`1e-5` (cf. :func:`~gudhi.DelaunayComplex.get_float_relative_precision`). - For more details, please refer to - `CGAL::Lazy_exact_nt::set_relative_precision_of_to_double `_ - :type precision: float + Args: + precision: When the DelaunayComplex is constructed with :code:`precision = 'safe'` (the default), + one can set the float relative precision of filtration values computed in + :func:`~gudhi.DelaunayComplex.create_simplex_tree`. + Default is :code:`1e-5` (cf. :func:`~gudhi.DelaunayComplex.get_float_relative_precision`). + For more details, please refer to + `CGAL::Lazy_exact_nt::set_relative_precision_of_to_double `_ """ if precision <=0. or precision >= 1.: raise ValueError("Relative precision value must be strictly greater than 0 and strictly lower than 1") Delaunay_complex_interface.set_float_relative_precision(precision) - + @staticmethod - def get_float_relative_precision(): + def get_float_relative_precision() -> float: """ - :returns: The float relative precision of filtration values computation in + Returns: + The float relative precision of filtration values computation in :func:`~gudhi.DelaunayComplex.create_simplex_tree` when the DelaunayComplex is constructed with :code:`precision = 'safe'` (the default). - :rtype: float """ return Delaunay_complex_interface.get_float_relative_precision() @@ -176,19 +179,21 @@ cdef class AlphaComplex(DelaunayComplex): When DelaunayComplex is constructed with an infinite value of alpha, the complex is a Delaunay complex. """ - def create_simplex_tree(self, max_alpha_square = float('inf'), default_filtration_value = False): + def create_simplex_tree(self, max_alpha_square: float = float('inf'), + default_filtration_value: bool = False, + square_root_filtrations: bool = False): """ - :param max_alpha_square: The maximum alpha square threshold the simplices shall not exceed. Default is set to - infinity, and there is very little point using anything else since it does not save time. - :type max_alpha_square: float - :param default_filtration_value: [Deprecated] Default value is `False` (which means compute the filtration - values). Set this value to `True` if filtration values are not needed to be computed (will be set to - `NaN`), but please consider constructing a :class:`~gudhi.DelaunayComplex` instead. - :type default_filtration_value: bool - :returns: A simplex tree created from the Delaunay Triangulation. The vertex `k` corresponds to the k-th input + Args: + max_alpha_square: The maximum alpha square threshold the simplices shall not exceed. Default is set to + infinity, and there is very little point using anything else since it does not save time. + default_filtration_value: [Deprecated] Default value is `False` (which means compute the filtration + values). Set this value to `True` if filtration values are not needed to be computed (will be set to + `NaN`), but please consider constructing a :class:`~gudhi.DelaunayComplex` instead. + square_root_filtrations: Square root filtration values when True. Default is False. + Returns: + SimplexTree: A simplex tree created from the Delaunay Triangulation. The vertex `k` corresponds to the k-th input point. The vertices may not be numbered contiguously as some points may be discarded in the triangulation (duplicate points, weighted hidden point, ...). - :rtype: SimplexTree """ filtration = 'alpha' if default_filtration_value: @@ -196,16 +201,16 @@ cdef class AlphaComplex(DelaunayComplex): warnings.warn('''Since Gudhi 3.10, creating an AlphaComplex with default_filtration_value=True is deprecated. Please consider constructing a DelaunayComplex instead. ''', DeprecationWarning) - return super().create_simplex_tree(max_alpha_square, filtration) + return super().create_simplex_tree(max_alpha_square, filtration, square_root_filtrations) - def get_point(self, vertex): + def get_point(self, vertex: int) -> list[float]: """This function returns the point corresponding to a given vertex from the :class:`~gudhi.SimplexTree` (the same as the k-th input point, where `k=vertex`) - :param vertex: The vertex. - :type vertex: int - :rtype: list of float - :returns: the point. + Args: + vertex: The vertex. + Returns: + the point. :raises IndexError: In case the point has no associated vertex in the diagram (because of weights or because it is a duplicate). @@ -225,25 +230,25 @@ cdef class DelaunayCechComplex(DelaunayComplex): When DelaunayCechComplex is constructed with an infinite value of alpha, the complex is a Delaunay complex. """ - def __init__(self, points=[], precision='safe'): + def __init__(self, points: Iterable[Iterable[float]] = [], precision: Literal['fast', 'safe', 'exact'] = 'safe'): """DelaunayCechComplex constructor. - :param points: A list of points in d-Dimension. - :type points: Iterable[Iterable[float]] - - :param precision: Delaunay Čech complex precision can be 'fast', 'safe' or 'exact'. Default is 'safe'. - :type precision: string + Args: + points (Iterable[Iterable[float]]): A list of points in d-Dimension. + precision (str): Delaunay Čech complex precision can be 'fast', 'safe' or 'exact'. Default is 'safe'. """ super().__init__(points = points, weights = [], precision = precision) - def create_simplex_tree(self, max_alpha_square = float('inf')): + def create_simplex_tree(self, max_alpha_square: float = float('inf'), + square_root_filtrations: bool = False): """ - :param max_alpha_square: The maximum alpha square threshold the simplices shall not exceed. Default is set to - infinity, and there is very little point using anything else since it does not save time. - :type max_alpha_square: float - :returns: A simplex tree created from the Delaunay Triangulation. The vertex `k` corresponds to the k-th input + Args: + max_alpha_square: The maximum alpha square threshold the simplices shall not exceed. Default is set to + infinity, and there is very little point using anything else since it does not save time. + square_root_filtrations: Square root filtration values when True. Default is False. + Returns: + SimplexTree: A simplex tree created from the Delaunay Triangulation. The vertex `k` corresponds to the k-th input point. The vertices may not be numbered contiguously as some points may be discarded in the triangulation (duplicate points, weighted hidden point, ...). - :rtype: SimplexTree """ - return super().create_simplex_tree(max_alpha_square, 'cech') + return super().create_simplex_tree(max_alpha_square, 'cech', square_root_filtrations) diff --git a/src/python/include/Delaunay_complex_factory.h b/src/python/include/Delaunay_complex_factory.h index 5265f6b697..aab042b2ac 100644 --- a/src/python/include/Delaunay_complex_factory.h +++ b/src/python/include/Delaunay_complex_factory.h @@ -7,6 +7,7 @@ * Modification(s): * - 2024/03 Vincent Rouvreau: Renamed Alpha_complex_factory as Delaunay_complex_factory for DelaunayCechComplex. * Factorize create_complex + * - 2024/10 Vincent Rouvreau: Add square root filtration values interface * - YYYY/MM Author: Description of the modification */ @@ -44,7 +45,7 @@ struct Point_cgal_to_cython; // Specialized Unweighted Functor template struct Point_cgal_to_cython { - std::vector operator()(CgalPointType const& point) const + std::vector operator()(CgalPointType const& point) const { std::vector vd; vd.reserve(point.dimension()); @@ -57,7 +58,7 @@ struct Point_cgal_to_cython { // Specialized Weighted Functor template struct Point_cgal_to_cython { - std::vector operator()(CgalPointType const& weighted_point) const + std::vector operator()(CgalPointType const& weighted_point) const { const auto& point = weighted_point.point(); return Point_cgal_to_cython()(point); @@ -73,7 +74,7 @@ static CgalPointType pt_cython_to_cgal(std::vector const& vec) { template bool create_complex(Delaunay_complex& delaunay_complex, Simplex_tree_interface* simplex_tree, const Point_cloud& points, double max_alpha_square, - bool exact_version, Delaunay_filtration filtration) { + bool exact_version, Delaunay_filtration filtration, bool square_root_filtrations) { if (filtration == Delaunay_filtration::CECH) { if (Weighted) throw std::runtime_error("Weighted Delaunay-Cech complex is not available"); @@ -84,13 +85,20 @@ bool create_complex(Delaunay_complex& delaunay_complex, Simplex_tree_interface* true); if (result == true) { // Construct the Delaunay-Cech complex by assigning filtration values with MEB - Gudhi::cech_complex::assign_MEB_filtration(Kernel(), *simplex_tree, points); + if (square_root_filtrations) + Gudhi::cech_complex::assign_MEB_filtration(Kernel(), *simplex_tree, points); + else + Gudhi::cech_complex::assign_MEB_filtration(Kernel(), *simplex_tree, points); simplex_tree->prune_above_filtration(max_alpha_square); } return result; } else { - return delaunay_complex.create_complex(*simplex_tree, max_alpha_square, - exact_version, filtration == Delaunay_filtration::NONE); + if (square_root_filtrations) + return delaunay_complex.template create_complex(*simplex_tree, max_alpha_square, + exact_version, filtration == Delaunay_filtration::NONE); + else + return delaunay_complex.template create_complex(*simplex_tree, max_alpha_square, + exact_version, filtration == Delaunay_filtration::NONE); } } @@ -100,10 +108,10 @@ class Abstract_delaunay_complex { virtual std::vector get_point(int vh) = 0; virtual bool create_simplex_tree(Simplex_tree_interface* simplex_tree, double max_alpha_square, - Delaunay_filtration filtration) = 0; - + Delaunay_filtration filtration, bool square_root_filtrations) = 0; + virtual std::size_t num_vertices() const = 0; - + virtual ~Abstract_delaunay_complex() = default; }; @@ -138,10 +146,11 @@ class Delaunay_complex_t final : public Abstract_delaunay_complex { } virtual bool create_simplex_tree(Simplex_tree_interface* simplex_tree, double max_alpha_square, - Delaunay_filtration filtration) override { + Delaunay_filtration filtration, bool square_root_filtrations) override { return create_complex>(delaunay_complex_, simplex_tree, points_, max_alpha_square, - exact_version_, filtration); + exact_version_, filtration, + square_root_filtrations); } virtual std::size_t num_vertices() const override { diff --git a/src/python/include/Delaunay_complex_interface.h b/src/python/include/Delaunay_complex_interface.h index 568f97b341..f44022e660 100644 --- a/src/python/include/Delaunay_complex_interface.h +++ b/src/python/include/Delaunay_complex_interface.h @@ -6,6 +6,7 @@ * * Modification(s): * - 2024/03 Vincent Rouvreau: Renamed Alpha_complex_interface as Delaunay_complex_interface for DelaunayCechComplex. + * - 2024/10 Vincent Rouvreau: Add square root filtration values interface * - YYYY/MM Author: Description of the modification */ @@ -81,10 +82,10 @@ class Delaunay_complex_interface { } void create_simplex_tree(Simplex_tree_interface* simplex_tree, double max_alpha_square, - Delaunay_filtration filtration) { + Delaunay_filtration filtration, bool square_root_filtrations) { // Nothing to be done in case of an empty point set if (delaunay_ptr_->num_vertices() > 0) - delaunay_ptr_->create_simplex_tree(simplex_tree, max_alpha_square, filtration); + delaunay_ptr_->create_simplex_tree(simplex_tree, max_alpha_square, filtration, square_root_filtrations); } static void set_float_relative_precision(double precision) { diff --git a/src/python/test/test_delaunay_complex.py b/src/python/test/test_delaunay_complex.py index 4becba648b..0583ccb16e 100644 --- a/src/python/test/test_delaunay_complex.py +++ b/src/python/test/test_delaunay_complex.py @@ -108,18 +108,18 @@ def _safe_persistence_comparison(simplicial_complex, precision): signal = [math.sin(x) for x in time] delta = math.pi delayed = [math.sin(x + delta) for x in time] - + #construct embedding embedding1 = [[signal[i], -signal[i]] for i in range(len(time))] embedding2 = [[signal[i], delayed[i]] for i in range(len(time))] - + #build simplicial_complex and simplex tree cplx1 = simplicial_complex(points=embedding1, precision = precision) simplex_tree1 = cplx1.create_simplex_tree() - + cplx2 = simplicial_complex(points=embedding2, precision = precision) simplex_tree2 = cplx2.create_simplex_tree() - + diag1 = simplex_tree1.persistence() diag2 = simplex_tree2.persistence() @@ -192,3 +192,15 @@ def test_duplicated_2d_points_on_a_plane(): for simplicial_complex in [AlphaComplex, DelaunayComplex, DelaunayCechComplex]: for precision in ['fast', 'safe', 'exact']: _duplicated_2d_points_on_a_plane(simplicial_complex, precision) + +def test_square_root_filtrations(): + for filtration in ['alpha', 'cech', None]: + for precision in ['fast', 'safe', 'exact']: + dc = DelaunayComplex(points=[[1, 1], [7, 0], [4, 6], [9, 6], [0, 14], [2, 19], [9, 17]], + precision = precision) + stree = dc.create_simplex_tree(filtration='cech', square_root_filtrations=False) + stree_sqrt = dc.create_simplex_tree(filtration='cech', square_root_filtrations=True) + for simplex, filt in stree_sqrt.get_filtration(): + # np.testing.assert_almost_equal(float('nan'), float('nan')) is ok + # while float('nan') == float('nan') is False + np.testing.assert_almost_equal(filt, math.sqrt(stree.filtration(simplex)))