From 986ad6f85d4d319ba957eca7d161e24f2fd3be6a Mon Sep 17 00:00:00 2001 From: Vincent Rouvreau Date: Tue, 1 Oct 2024 12:34:15 +0200 Subject: [PATCH 01/17] c++ part of square root filtration values for alpha and delaunay-cech complex --- .../include/gudhi/Alpha_complex.h | 50 +++++++++++-------- .../test/Delaunay_complex_unit_test.h | 18 ++++++- .../test/Weighted_alpha_complex_unit_test.cpp | 45 ++++++++++++++++- src/Alpha_complex/utilities/CMakeLists.txt | 20 +++++++- .../utilities/alpha_complex_persistence.cpp | 50 ++++++++++++------- src/Alpha_complex/utilities/alphacomplex.md | 14 ++++-- .../include/gudhi/MEB_filtration.h | 38 ++++++++++---- src/Cech_complex/test/test_cech_complex.cpp | 24 +++++++++ src/Cech_complex/utilities/CMakeLists.txt | 18 +++++++ src/Cech_complex/utilities/cechcomplex.md | 10 +++- .../utilities/delaunay_cech_persistence.cpp | 30 ++++++++--- 11 files changed, 253 insertions(+), 64 deletions(-) diff --git a/src/Alpha_complex/include/gudhi/Alpha_complex.h b/src/Alpha_complex/include/gudhi/Alpha_complex.h index ca9962046d..24e5dd055f 100644 --- a/src/Alpha_complex/include/gudhi/Alpha_complex.h +++ b/src/Alpha_complex/include/gudhi/Alpha_complex.h @@ -65,14 +65,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 +161,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 +183,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 +198,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. */ @@ -355,7 +355,7 @@ class Alpha_complex { * 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. Useless if `default_filtration_value` is set to @@ -367,13 +367,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 +390,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 +442,6 @@ class Alpha_complex { } } // -------------------------------------------------------------------------------------------- - if (!default_filtration_value) { CGAL::NT_converter cgal_converter; // -------------------------------------------------------------------------------------------- @@ -458,6 +461,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 +479,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..f6583f9328 100644 --- a/src/Alpha_complex/utilities/alpha_complex_persistence.cpp +++ b/src/Alpha_complex/utilities/alpha_complex_persistence.cpp @@ -29,8 +29,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 +62,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 +78,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 +104,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 +122,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 +155,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 +169,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 +202,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..50dc86e46b 100644 --- a/src/Cech_complex/utilities/delaunay_cech_persistence.cpp +++ b/src/Cech_complex/utilities/delaunay_cech_persistence.cpp @@ -31,15 +31,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 +51,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 +66,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 +84,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 +117,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 +130,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 +161,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"; From 6be76bc4e837dd33d620386c56df2cccd75befec Mon Sep 17 00:00:00 2001 From: Vincent Rouvreau Date: Tue, 1 Oct 2024 14:31:14 +0200 Subject: [PATCH 02/17] square_root_filtrations was missing in documentation for Alpha complex --- src/Alpha_complex/include/gudhi/Alpha_complex.h | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/Alpha_complex/include/gudhi/Alpha_complex.h b/src/Alpha_complex/include/gudhi/Alpha_complex.h index 24e5dd055f..bffb6a1c53 100644 --- a/src/Alpha_complex/include/gudhi/Alpha_complex.h +++ b/src/Alpha_complex/include/gudhi/Alpha_complex.h @@ -354,6 +354,8 @@ 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. From 7b14133413699873165139854759d004d1a04b8e Mon Sep 17 00:00:00 2001 From: Vincent Rouvreau Date: Wed, 2 Oct 2024 11:35:36 +0200 Subject: [PATCH 03/17] Add new feature in changes --- src/Alpha_complex/include/gudhi/Alpha_complex.h | 1 + src/Alpha_complex/utilities/alpha_complex_persistence.cpp | 1 + src/Cech_complex/utilities/delaunay_cech_persistence.cpp | 1 + 3 files changed, 3 insertions(+) diff --git a/src/Alpha_complex/include/gudhi/Alpha_complex.h b/src/Alpha_complex/include/gudhi/Alpha_complex.h index bffb6a1c53..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 */ diff --git a/src/Alpha_complex/utilities/alpha_complex_persistence.cpp b/src/Alpha_complex/utilities/alpha_complex_persistence.cpp index f6583f9328..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 */ diff --git a/src/Cech_complex/utilities/delaunay_cech_persistence.cpp b/src/Cech_complex/utilities/delaunay_cech_persistence.cpp index 50dc86e46b..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 */ From 2cdc152899b5b2795ba12b754e7d107d77873c6d Mon Sep 17 00:00:00 2001 From: Vincent Rouvreau Date: Wed, 2 Oct 2024 15:40:10 +0200 Subject: [PATCH 04/17] Python version of square root filtrations for alpha complex and delaunay cech. Python hints for delaunay complex family --- src/python/gudhi/delaunay_complex.pyx | 134 +++++++++--------- src/python/include/Delaunay_complex_factory.h | 31 ++-- .../include/Delaunay_complex_interface.h | 5 +- 3 files changed, 92 insertions(+), 78 deletions(-) diff --git a/src/python/gudhi/delaunay_complex.pyx b/src/python/gudhi/delaunay_complex.pyx index 49144c89fa..713cbc017e 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 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: A list of points in d-Dimension. + weights: A list of weights. If set, the number of weights must correspond to the number of points. + precision: 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,30 @@ 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 +178,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 +200,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 +229,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: A list of points in d-Dimension. + precision: 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) { From c578b93baee36462a0d540128d975f7b5edd5319 Mon Sep 17 00:00:00 2001 From: Vincent Rouvreau Date: Thu, 3 Oct 2024 10:30:14 +0200 Subject: [PATCH 05/17] python hints not well managed on cython class --- src/python/gudhi/delaunay_complex.pyx | 15 ++++++++------- 1 file changed, 8 insertions(+), 7 deletions(-) diff --git a/src/python/gudhi/delaunay_complex.pyx b/src/python/gudhi/delaunay_complex.pyx index 713cbc017e..8938b57e6b 100644 --- a/src/python/gudhi/delaunay_complex.pyx +++ b/src/python/gudhi/delaunay_complex.pyx @@ -19,7 +19,7 @@ from libcpp.string cimport string from libcpp cimport bool from libc.stdint cimport intptr_t import warnings -from typing import Literal, Optional +from typing import Literal, Optional, Iterable from gudhi.simplex_tree cimport * from gudhi.simplex_tree import SimplexTree @@ -67,9 +67,9 @@ cdef class DelaunayComplex: """DelaunayComplex constructor. Args: - points: A list of points in d-Dimension. - weights: A list of weights. If set, the number of weights must correspond to the number of points. - precision: Delaunay complex precision can be 'fast', 'safe' or 'exact'. Default is 'safe'. + 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. """ @@ -155,7 +155,8 @@ cdef class DelaunayComplex: @staticmethod 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). """ @@ -233,8 +234,8 @@ cdef class DelaunayCechComplex(DelaunayComplex): """DelaunayCechComplex constructor. Args: - points: A list of points in d-Dimension. - precision: Delaunay Čech complex precision can be 'fast', 'safe' or 'exact'. Default is 'safe'. + 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) From 46f5616debc3ca02de82e54fe52ea90d5db64145 Mon Sep 17 00:00:00 2001 From: Vincent Rouvreau Date: Mon, 14 Oct 2024 13:15:04 +0200 Subject: [PATCH 06/17] Add some python test for square_root_filtrations versions --- src/python/test/test_delaunay_complex.py | 20 ++++++++++++++++---- 1 file changed, 16 insertions(+), 4 deletions(-) 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))) From 4d14cbcf850532c8114c75e9c3c6334f3cfe8f4a Mon Sep 17 00:00:00 2001 From: Vincent Rouvreau Date: Thu, 21 Nov 2024 18:43:15 +0100 Subject: [PATCH 07/17] code review: Test filtration values are positive, before sqrt --- src/Cech_complex/include/gudhi/MEB_filtration.h | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/Cech_complex/include/gudhi/MEB_filtration.h b/src/Cech_complex/include/gudhi/MEB_filtration.h index cd30a7ead4..f0c2ed7505 100644 --- a/src/Cech_complex/include/gudhi/MEB_filtration.h +++ b/src/Cech_complex/include/gudhi/MEB_filtration.h @@ -81,10 +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()); - Filtration_value filt{cvt(r)}; + Filtration_value filt{max(cvt(r), Filtration_value(0))}; if constexpr (square_root_filtrations) filt = sqrt(filt); - complex.assign_filtration(sh, max(filt, Filtration_value(0))); + complex.assign_filtration(sh, filt); 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 @@ -121,7 +121,7 @@ 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. - Filtration_value filt{cvt(r)}; + Filtration_value filt{max(cvt(r), Filtration_value(0))}; if constexpr (square_root_filtrations) filt = sqrt(filt); // maxf = filt except for rounding errors From f80347f8812537db879a2c2cbb0bc8a18341601d Mon Sep 17 00:00:00 2001 From: Vincent Rouvreau <10407034+VincentRouvreau@users.noreply.github.com> Date: Thu, 21 Nov 2024 18:52:11 +0100 Subject: [PATCH 08/17] code review: typing.Iterable s deprecated since python 3.9 Co-authored-by: David Loiseaux --- src/python/gudhi/delaunay_complex.pyx | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/python/gudhi/delaunay_complex.pyx b/src/python/gudhi/delaunay_complex.pyx index 8938b57e6b..31662b9567 100644 --- a/src/python/gudhi/delaunay_complex.pyx +++ b/src/python/gudhi/delaunay_complex.pyx @@ -19,7 +19,8 @@ from libcpp.string cimport string from libcpp cimport bool from libc.stdint cimport intptr_t import warnings -from typing import Literal, Optional, Iterable +from typing import Literal, Optional +from collections.abc import Iterable from gudhi.simplex_tree cimport * from gudhi.simplex_tree import SimplexTree From dc78dc15bb324fcb0c8c9f4eb43115755cb11589 Mon Sep 17 00:00:00 2001 From: Vincent Rouvreau Date: Fri, 22 Nov 2024 11:33:47 +0100 Subject: [PATCH 09/17] mypy + cython + sphinx compromise --- src/python/gudhi/delaunay_complex.pyx | 16 +++++++--------- 1 file changed, 7 insertions(+), 9 deletions(-) diff --git a/src/python/gudhi/delaunay_complex.pyx b/src/python/gudhi/delaunay_complex.pyx index 31662b9567..c3891e0b6f 100644 --- a/src/python/gudhi/delaunay_complex.pyx +++ b/src/python/gudhi/delaunay_complex.pyx @@ -104,9 +104,9 @@ cdef class DelaunayComplex: """ return self.this_ptr != NULL - def create_simplex_tree(self, max_alpha_square: float = float('inf'), + def create_simplex_tree(self, double max_alpha_square: float = float('inf'), filtration: Optional[Literal['alpha', 'cech']] = None, - square_root_filtrations: bool = False): + bool square_root_filtrations: bool = False) -> SimplexTree: """ Args: max_alpha_square: The maximum alpha square threshold the simplices shall not exceed. Default is set to @@ -123,9 +123,7 @@ cdef class DelaunayComplex: 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': @@ -135,7 +133,7 @@ cdef class DelaunayComplex: with nogil: self.this_ptr.create_simplex_tree(stree_int_ptr, - mas, filt, srf) + max_alpha_square, filt, square_root_filtrations) return stree @staticmethod @@ -180,9 +178,9 @@ 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 = float('inf'), + def create_simplex_tree(self, double max_alpha_square: float = float('inf'), default_filtration_value: bool = False, - square_root_filtrations: bool = False): + bool square_root_filtrations: bool = False) -> SimplexTree: """ Args: max_alpha_square: The maximum alpha square threshold the simplices shall not exceed. Default is set to @@ -240,8 +238,8 @@ cdef class DelaunayCechComplex(DelaunayComplex): """ super().__init__(points = points, weights = [], precision = precision) - def create_simplex_tree(self, max_alpha_square: float = float('inf'), - square_root_filtrations: bool = False): + def create_simplex_tree(self, double max_alpha_square: float = float('inf'), + bool square_root_filtrations: bool = False) -> SimplexTree: """ Args: max_alpha_square: The maximum alpha square threshold the simplices shall not exceed. Default is set to From 5d9a56463223d199d82e1fb15b08449348a368fa Mon Sep 17 00:00:00 2001 From: Vincent Rouvreau Date: Fri, 22 Nov 2024 12:09:19 +0100 Subject: [PATCH 10/17] code review: rename square_root_filtrations as output_squared_values --- .../include/gudhi/Alpha_complex.h | 10 +++++----- .../test/Delaunay_complex_unit_test.h | 2 +- .../test/Weighted_alpha_complex_unit_test.cpp | 2 +- .../utilities/alpha_complex_persistence.cpp | 20 +++++++++---------- .../include/gudhi/MEB_filtration.h | 16 +++++++-------- src/Cech_complex/test/test_cech_complex.cpp | 4 ++-- .../utilities/delaunay_cech_persistence.cpp | 20 +++++++++---------- src/python/gudhi/delaunay_complex.pyx | 20 +++++++++---------- src/python/include/Delaunay_complex_factory.h | 12 +++++------ .../include/Delaunay_complex_interface.h | 4 ++-- src/python/test/test_delaunay_complex.py | 6 +++--- 11 files changed, 58 insertions(+), 58 deletions(-) diff --git a/src/Alpha_complex/include/gudhi/Alpha_complex.h b/src/Alpha_complex/include/gudhi/Alpha_complex.h index 3074858db1..e9b9ab8c43 100644 --- a/src/Alpha_complex/include/gudhi/Alpha_complex.h +++ b/src/Alpha_complex/include/gudhi/Alpha_complex.h @@ -355,8 +355,8 @@ 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 output_squared_values 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 output_squared_values is true. * \tparam SimplicialComplexForAlpha must meet `SimplicialComplexForAlpha` concept. * * @param[in] complex SimplicialComplexForAlpha to be created. @@ -376,7 +376,7 @@ class Alpha_complex { * * Initialization can be launched once. */ - template bool create_complex(SimplicialComplexForAlpha& complex, @@ -464,7 +464,7 @@ class Alpha_complex { if(exact) CGAL::exact(sqrad); #endif alpha_complex_filtration = cgal_converter(sqrad); - if constexpr (square_root_filtrations) { + if constexpr (output_squared_values) { alpha_complex_filtration = sqrt(alpha_complex_filtration); } } @@ -489,7 +489,7 @@ class Alpha_complex { // 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 - if constexpr (square_root_filtrations) { + if constexpr (output_squared_values) { complex.prune_above_filtration(sqrt(max_alpha_square)); } else { complex.prune_above_filtration(max_alpha_square); diff --git a/src/Alpha_complex/test/Delaunay_complex_unit_test.h b/src/Alpha_complex/test/Delaunay_complex_unit_test.h index df24441f6e..aa16423b1a 100644 --- a/src/Alpha_complex/test/Delaunay_complex_unit_test.h +++ b/src/Alpha_complex/test/Delaunay_complex_unit_test.h @@ -51,7 +51,7 @@ void compare_delaunay_complex_simplices() { 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"; + std::clog << "Check simplices from alpha complex filtration values when output_squared_values 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)); 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 f40521761a..8a353228fc 100644 --- a/src/Alpha_complex/test/Weighted_alpha_complex_unit_test.cpp +++ b/src/Alpha_complex/test/Weighted_alpha_complex_unit_test.cpp @@ -154,7 +154,7 @@ BOOST_AUTO_TEST_CASE(Is_weighted_alpha_complex_nan) { BOOST_CHECK(!std::isnan(stree.filtration(f_simplex))); } } - std::clog << "Weighted alpha complex with square_root_filtrations\n"; + std::clog << "Weighted alpha complex with output_squared_values\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()) { diff --git a/src/Alpha_complex/utilities/alpha_complex_persistence.cpp b/src/Alpha_complex/utilities/alpha_complex_persistence.cpp index 5cc4e997c5..b8cc6fc19a 100644 --- a/src/Alpha_complex/utilities/alpha_complex_persistence.cpp +++ b/src/Alpha_complex/utilities/alpha_complex_persistence.cpp @@ -30,7 +30,7 @@ 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, - bool& square_root_filtrations, std::string &weight_file, std::string &output_file_diag, + bool& output_squared_values, std::string &weight_file, std::string &output_file_diag, Filtration_value &alpha_square_max_value, int &coeff_field_characteristic, Filtration_value &min_persistence); @@ -64,7 +64,7 @@ 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 square_root_filtrations) { + bool output_squared_values) { Simplex_tree stree; auto points = read_off(off_file_points); @@ -79,7 +79,7 @@ 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 (square_root_filtrations) + if (output_squared_values) 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); @@ -87,7 +87,7 @@ 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); - if (square_root_filtrations) + if (output_squared_values) 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); @@ -105,12 +105,12 @@ int main(int argc, char **argv) { std::string output_file_diag; bool exact_version = false; bool fast_version = false; - bool square_root_filtrations = false; + bool output_squared_values = 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, square_root_filtrations, weight_file, + program_options(argc, argv, off_file_points, exact_version, fast_version, output_squared_values, weight_file, output_file_diag, alpha_square_max_value, coeff_field_characteristic, min_persistence); if ((exact_version) && (fast_version)) { @@ -123,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, square_root_filtrations); + stree = create_simplex_tree(off_file_points, weight_file, exact_version, alpha_square_max_value, output_squared_values); } else { using Kernel = CGAL::Epeck_d; - stree = create_simplex_tree(off_file_points, weight_file, exact_version, alpha_square_max_value, square_root_filtrations); + stree = create_simplex_tree(off_file_points, weight_file, exact_version, alpha_square_max_value, output_squared_values); } // ---------------------------------------------------------------------------- // Display information about the alpha complex @@ -156,7 +156,7 @@ int main(int argc, char **argv) { } void program_options(int argc, char *argv[], std::string &off_file_points, bool &exact, bool &fast, - bool& square_root_filtrations, std::string &weight_file, std::string &output_file_diag, + bool& output_squared_values, 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; @@ -170,7 +170,7 @@ 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), + "square-root-filtrations,s", po::bool_switch(&output_squared_values), "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 ")( diff --git a/src/Cech_complex/include/gudhi/MEB_filtration.h b/src/Cech_complex/include/gudhi/MEB_filtration.h index f0c2ed7505..6f61f68664 100644 --- a/src/Cech_complex/include/gudhi/MEB_filtration.h +++ b/src/Cech_complex/include/gudhi/MEB_filtration.h @@ -5,7 +5,7 @@ * Copyright (C) 2023 Inria * * Modification(s): - * - 2024/10 Vincent Rouvreau: Add square_root_filtrations argument to enable/disable squared radii computation + * - 2024/10 Vincent Rouvreau: Add output_squared_values argument to enable/disable squared radii computation * - YYYY/MM Author: Description of the modification */ @@ -25,13 +25,13 @@ namespace Gudhi::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 (or not squared in function of `square_root_filtrations`) radius of its minimal enclosing ball (MEB). + * to the squared (or not squared in function of `output_squared_values`) radius of its minimal enclosing ball (MEB). * - * Applied on a Čech complex, it recomputes the same values (squared or not in function of `square_root_filtrations`). + * Applied on a Čech complex, it recomputes the same values (squared or not in function of `output_squared_values`). * 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 output_squared_values 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 output_squared_values is true. * \tparam Kernel CGAL kernel: either Epick_d or Epeck_d. * \tparam PointRange Random access range of `Kernel::Point_d`. * @@ -42,7 +42,7 @@ namespace Gudhi::cech_complex { * 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; @@ -82,7 +82,7 @@ void assign_MEB_filtration(Kernel&&k, SimplicialComplexForMEB& complex, PointRan if (exact) CGAL::exact(r); complex.assign_key(sh, cache_.size()); Filtration_value filt{max(cvt(r), Filtration_value(0))}; - if constexpr (square_root_filtrations) + if constexpr (output_squared_values) filt = sqrt(filt); complex.assign_filtration(sh, filt); cache_.emplace_back(std::move(m), std::move(r)); @@ -122,7 +122,7 @@ void assign_MEB_filtration(Kernel&&k, SimplicialComplexForMEB& complex, PointRan // Filtration_value max_sanity = maxf * d2 / (d2 - 1); // and use min(max_sanity, ...), which would limit how bad numerical errors can be. Filtration_value filt{max(cvt(r), Filtration_value(0))}; - if constexpr (square_root_filtrations) + if constexpr (output_squared_values) filt = sqrt(filt); // maxf = filt except for rounding errors maxf = max(maxf, filt); diff --git a/src/Cech_complex/test/test_cech_complex.cpp b/src/Cech_complex/test/test_cech_complex.cpp index 8d0901834e..061355d6ce 100644 --- a/src/Cech_complex/test/test_cech_complex.cpp +++ b/src/Cech_complex/test/test_cech_complex.cpp @@ -193,10 +193,10 @@ BOOST_AUTO_TEST_CASE(Cech_complex_for_documentation) { } BOOST_CHECK(st3 == st3_save); // Should only be an approximate test - // Same test but with square_root_filtrations=true + // Same test but with output_squared_values=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"; + std::clog << "Simplex tree from assign_MEB_filtration with output_squared_values=true\n"; for (auto f_simplex : st3.filtration_simplex_range()) { std::clog << " ( "; for (auto vertex : st3.simplex_vertex_range(f_simplex)) std::clog << vertex << " "; diff --git a/src/Cech_complex/utilities/delaunay_cech_persistence.cpp b/src/Cech_complex/utilities/delaunay_cech_persistence.cpp index e967a2b5bb..fa20854dea 100644 --- a/src/Cech_complex/utilities/delaunay_cech_persistence.cpp +++ b/src/Cech_complex/utilities/delaunay_cech_persistence.cpp @@ -5,7 +5,7 @@ * Copyright (C) 2024 Inria * * Modification(s): - * - 2024/10 Vincent Rouvreau: Add square_root_filtrations argument to enable/disable squared radii computation + * - 2024/10 Vincent Rouvreau: Add output_squared_values argument to enable/disable squared radii computation * - YYYY/MM Author: Description of the modification */ @@ -32,12 +32,12 @@ 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, - bool& square_root_filtrations, std::string& pers_diag_file, Filtration_value& max_radius, int& p, + bool& output_squared_values, 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, - bool square_root_filtrations) { + bool output_squared_values) { using Point = typename Kernel::Point_d; using Points_off_reader = Gudhi::Points_off_reader; using Delaunay_complex = Gudhi::alpha_complex::Alpha_complex; @@ -52,7 +52,7 @@ 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); - if (square_root_filtrations) { + if (output_squared_values) { Gudhi::cech_complex::assign_MEB_filtration(Kernel(), stree, point_cloud, exact_version); max_radius = sqrt(max_radius); } else { @@ -67,12 +67,12 @@ int main(int argc, char* argv[]) { std::string pers_diag_file; bool exact_version = false; bool fast_version = false; - bool square_root_filtrations = false; + bool output_squared_values = false; Filtration_value max_radius; int p; Filtration_value min_persistence; - program_options(argc, argv, off_file_points, exact_version, fast_version, square_root_filtrations, pers_diag_file, + program_options(argc, argv, off_file_points, exact_version, fast_version, output_squared_values, pers_diag_file, max_radius, p, min_persistence); if ((exact_version) && (fast_version)) { @@ -85,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, square_root_filtrations); + stree = create_simplex_tree(off_file_points, exact_version, max_radius, output_squared_values); } else { std::clog << "exact_version = " << exact_version << "\n"; using Kernel = CGAL::Epeck_d; - stree = create_simplex_tree(off_file_points, exact_version, max_radius, square_root_filtrations); + stree = create_simplex_tree(off_file_points, exact_version, max_radius, output_squared_values); } std::clog << "The complex contains " << stree.num_simplices() << " simplices \n"; @@ -118,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, - bool& square_root_filtrations, std::string& pers_diag_file, Filtration_value& max_radius, int& p, + bool& output_squared_values, 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"); @@ -131,7 +131,7 @@ 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), + "square-root-filtrations,s", po::bool_switch(&output_squared_values), "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")( diff --git a/src/python/gudhi/delaunay_complex.pyx b/src/python/gudhi/delaunay_complex.pyx index c3891e0b6f..8989c75f0b 100644 --- a/src/python/gudhi/delaunay_complex.pyx +++ b/src/python/gudhi/delaunay_complex.pyx @@ -43,7 +43,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, bool square_root_filtrations) nogil except + + void create_simplex_tree(Simplex_tree_python_interface* simplex_tree, double max_alpha_square, Delaunay_filtration filtration, bool output_squared_values) nogil except + @staticmethod void set_float_relative_precision(double precision) nogil @staticmethod @@ -106,7 +106,7 @@ cdef class DelaunayComplex: def create_simplex_tree(self, double max_alpha_square: float = float('inf'), filtration: Optional[Literal['alpha', 'cech']] = None, - bool square_root_filtrations: bool = False) -> SimplexTree: + bool output_squared_values: bool = False) -> SimplexTree: """ Args: max_alpha_square: The maximum alpha square threshold the simplices shall not exceed. Default is set to @@ -114,7 +114,7 @@ cdef class DelaunayComplex: 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. + output_squared_values: 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 @@ -133,7 +133,7 @@ cdef class DelaunayComplex: with nogil: self.this_ptr.create_simplex_tree(stree_int_ptr, - max_alpha_square, filt, square_root_filtrations) + max_alpha_square, filt, output_squared_values) return stree @staticmethod @@ -180,7 +180,7 @@ cdef class AlphaComplex(DelaunayComplex): """ def create_simplex_tree(self, double max_alpha_square: float = float('inf'), default_filtration_value: bool = False, - bool square_root_filtrations: bool = False) -> SimplexTree: + bool output_squared_values: bool = False) -> SimplexTree: """ Args: max_alpha_square: The maximum alpha square threshold the simplices shall not exceed. Default is set to @@ -188,7 +188,7 @@ cdef class AlphaComplex(DelaunayComplex): 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. + output_squared_values: 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 @@ -200,7 +200,7 @@ 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, square_root_filtrations) + return super().create_simplex_tree(max_alpha_square, filtration, output_squared_values) def get_point(self, vertex: int) -> list[float]: """This function returns the point corresponding to a given vertex from the :class:`~gudhi.SimplexTree` (the @@ -239,15 +239,15 @@ cdef class DelaunayCechComplex(DelaunayComplex): super().__init__(points = points, weights = [], precision = precision) def create_simplex_tree(self, double max_alpha_square: float = float('inf'), - bool square_root_filtrations: bool = False) -> SimplexTree: + bool output_squared_values: bool = False) -> 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. - square_root_filtrations: Square root filtration values when True. Default is False. + output_squared_values: 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, ...). """ - return super().create_simplex_tree(max_alpha_square, 'cech', square_root_filtrations) + return super().create_simplex_tree(max_alpha_square, 'cech', output_squared_values) diff --git a/src/python/include/Delaunay_complex_factory.h b/src/python/include/Delaunay_complex_factory.h index aab042b2ac..c945d1637d 100644 --- a/src/python/include/Delaunay_complex_factory.h +++ b/src/python/include/Delaunay_complex_factory.h @@ -74,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 square_root_filtrations) { + bool exact_version, Delaunay_filtration filtration, bool output_squared_values) { if (filtration == Delaunay_filtration::CECH) { if (Weighted) throw std::runtime_error("Weighted Delaunay-Cech complex is not available"); @@ -85,7 +85,7 @@ 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 - if (square_root_filtrations) + if (output_squared_values) Gudhi::cech_complex::assign_MEB_filtration(Kernel(), *simplex_tree, points); else Gudhi::cech_complex::assign_MEB_filtration(Kernel(), *simplex_tree, points); @@ -93,7 +93,7 @@ bool create_complex(Delaunay_complex& delaunay_complex, Simplex_tree_interface* } return result; } else { - if (square_root_filtrations) + if (output_squared_values) return delaunay_complex.template create_complex(*simplex_tree, max_alpha_square, exact_version, filtration == Delaunay_filtration::NONE); else @@ -108,7 +108,7 @@ 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, bool square_root_filtrations) = 0; + Delaunay_filtration filtration, bool output_squared_values) = 0; virtual std::size_t num_vertices() const = 0; @@ -146,11 +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, bool square_root_filtrations) override { + Delaunay_filtration filtration, bool output_squared_values) override { return create_complex>(delaunay_complex_, simplex_tree, points_, max_alpha_square, exact_version_, filtration, - square_root_filtrations); + output_squared_values); } 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 f44022e660..5fdb8970e4 100644 --- a/src/python/include/Delaunay_complex_interface.h +++ b/src/python/include/Delaunay_complex_interface.h @@ -82,10 +82,10 @@ class Delaunay_complex_interface { } void create_simplex_tree(Simplex_tree_interface* simplex_tree, double max_alpha_square, - Delaunay_filtration filtration, bool square_root_filtrations) { + Delaunay_filtration filtration, bool output_squared_values) { // 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, square_root_filtrations); + delaunay_ptr_->create_simplex_tree(simplex_tree, max_alpha_square, filtration, output_squared_values); } 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 0583ccb16e..84281fcd11 100644 --- a/src/python/test/test_delaunay_complex.py +++ b/src/python/test/test_delaunay_complex.py @@ -193,13 +193,13 @@ def test_duplicated_2d_points_on_a_plane(): for precision in ['fast', 'safe', 'exact']: _duplicated_2d_points_on_a_plane(simplicial_complex, precision) -def test_square_root_filtrations(): +def test_output_squared_values(): 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) + stree = dc.create_simplex_tree(filtration='cech', output_squared_values=False) + stree_sqrt = dc.create_simplex_tree(filtration='cech', output_squared_values=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 From f3ca164a5a93eedf372d79c5ac661e053e435da3 Mon Sep 17 00:00:00 2001 From: Vincent Rouvreau Date: Tue, 26 Nov 2024 11:43:12 +0100 Subject: [PATCH 11/17] code review: go with free function instead of classes --- src/python/doc/delaunay_complex_ref.rst | 34 ++- src/python/doc/delaunay_complex_user.rst | 102 ++++---- src/python/gudhi/delaunay_complex.pyx | 187 +++++++++++---- src/python/include/Delaunay_complex_factory.h | 9 +- src/python/test/test_delaunay_complex.py | 218 ++++++++++-------- 5 files changed, 338 insertions(+), 212 deletions(-) diff --git a/src/python/doc/delaunay_complex_ref.rst b/src/python/doc/delaunay_complex_ref.rst index b4e35a80d6..fb18ae047a 100644 --- a/src/python/doc/delaunay_complex_ref.rst +++ b/src/python/doc/delaunay_complex_ref.rst @@ -6,21 +6,35 @@ Delaunay complex reference manual ================================= -.. autoclass:: gudhi.DelaunayComplex - :members: +.. autofunction:: gudhi.delaunay_complex + +====================================== +Delaunay Čech complex reference manual +====================================== + +.. autofunction:: gudhi.delaunay_cech_complex ============================== Alpha complex reference manual ============================== +.. autofunction:: gudhi.alpha_complex + +.. autofunction:: gudhi.weighted_alpha_complex + .. autoclass:: gudhi.AlphaComplex - :members: - :show-inheritance: -====================================== -Delaunay Čech complex reference manual -====================================== + .. deprecated:: 3.11.0 + Use :py:func:`gudhi.alpha_complex`, :py:func:`gudhi.weighted_alpha_complex` or :py:func:`gudhi.delaunay_complex` + instead. + + .. automethod:: gudhi.AlphaComplex.create_simplex_tree + .. automethod:: gudhi.AlphaComplex.get_point + +================== +Relative precision +================== + +.. autofunction:: gudhi.DelaunayComplex.set_float_relative_precision -.. autoclass:: gudhi.DelaunayCechComplex - :members: - :show-inheritance: \ No newline at end of file +.. autofunction:: gudhi.DelaunayComplex.get_float_relative_precision diff --git a/src/python/doc/delaunay_complex_user.rst b/src/python/doc/delaunay_complex_user.rst index 844ccd609c..fb2e966c67 100644 --- a/src/python/doc/delaunay_complex_user.rst +++ b/src/python/doc/delaunay_complex_user.rst @@ -10,26 +10,21 @@ Definition .. include:: delaunay_complex_sum.inc -:class:`~gudhi.DelaunayComplex` is constructing a :doc:`SimplexTree ` using +:func:`~gudhi.delaunay_complex` is constructing a :doc:`SimplexTree ` using `Delaunay Triangulation `_ :cite:`cgal:hdj-t-19b` from the `Computational Geometry Algorithms Library `_ -:cite:`cgal:eb-19b`. +:cite:`cgal:eb-19b`. In this case, all filtration values are set to `NaN` -The Delaunay complex (all filtration values are set to `NaN`) is available by passing :code:`filtrations = None` -(default value) to :func:`~gudhi.DelaunayComplex.create_simplex_tree`. +In order to compute filtration values, :func:`~gudhi.delaunay_cech_complex` is a bit faster than +:func:`~gudhi.alpha_complex`, and a weighted version, :func:`~gudhi.weighted_alpha_complex` is available. -When :paramref:`~gudhi.DelaunayComplex.create_simplex_tree.filtrations` is: +Remarks about the filtration value computation +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ -* `'alpha'` - The filtration value of each simplex is computed as an :class:`~gudhi.AlphaComplex` -* `'cech'` - The filtration value of each simplex is computed as a :class:`~gudhi.DelaunayCechComplex` +These remarks apply to the :func:`~gudhi.delaunay_cech_complex`, :func:`~gudhi.alpha_complex`, and +:func:`~gudhi.weighted_alpha_complex` that we will call *'complex'* in the following text. -Remarks about the :class:`~gudhi.AlphaComplex` and the :class:`~gudhi.DelaunayCechComplex` -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -These remarks apply to the :class:`~gudhi.AlphaComplex` and to the :class:`~gudhi.DelaunayCechComplex` that we will -call *'complex'* in the following text. - -* When a complex is constructed with an infinite value of :math:`\alpha^2`, the complex is a Delaunay complex (with +* When a complex is constructed with an infinite value of :math:`\alpha`, the complex is a Delaunay complex (with special filtration values). * For people only interested in the topology of the complex (for instance persistence), notice the complex is equivalent to the `Čech complex `_ and much smaller if @@ -46,27 +41,22 @@ call *'complex'* in the following text. the filtration values can exceptionally be arbitrarily bad. In all cases, we still guarantee that the output is a valid filtration (faces have a filtration value no larger than their cofaces). -:class:`~gudhi.DelaunayCechComplex` is a bit faster than :class:`~gudhi.AlphaComplex`, but only -:class:`~gudhi.AlphaComplex` has a weighted version. - Example from points ------------------- -This example builds the Delaunay Čech complex from the given points: +This example builds the Delaunay Čech complex (represented as a :doc:`SimplexTree `) from the given +points: .. testcode:: - from gudhi import DelaunayCechComplex + from gudhi import delaunay_cech_complex points=[[1, 1], [7, 0], [4, 6], [9, 6], [0, 14], [2, 19], [9, 17]] - cpx = DelaunayCechComplex(points=points) + stree = delaunay_cech_complex(points=points) - stree = cpx.create_simplex_tree() - print('Complex is of dimension ', stree.dimension(), ' - ', - stree.num_simplices(), ' simplices - ', stree.num_vertices(), ' vertices.') - - fmt = '%s -> %.2f' + print(f"Complex is of dimension {stree.dimension()} - {stree.num_simplices()} simplices - ", + f"{stree.num_vertices()} vertices.") for filtered_value in stree.get_filtration(): - print(fmt % tuple(filtered_value)) + print('%s -> %.2f' % tuple(filtered_value)) The output is: @@ -81,24 +71,24 @@ The output is: [4] -> 0.00 [5] -> 0.00 [6] -> 0.00 - [2, 3] -> 6.25 - [4, 5] -> 7.25 - [0, 2] -> 8.50 - [0, 1] -> 9.25 - [1, 3] -> 10.00 - [1, 2] -> 11.25 - [1, 2, 3] -> 12.50 - [0, 1, 2] -> 13.00 - [5, 6] -> 13.25 - [2, 4] -> 20.00 - [4, 6] -> 22.50 - [4, 5, 6] -> 22.50 - [3, 6] -> 30.25 - [2, 6] -> 36.50 - [2, 3, 6] -> 36.50 - [2, 4, 6] -> 37.24 - [0, 4] -> 42.50 - [0, 2, 4] -> 42.50 + [2, 3] -> 2.50 + [4, 5] -> 2.69 + [0, 2] -> 2.92 + [0, 1] -> 3.04 + [1, 3] -> 3.16 + [1, 2] -> 3.35 + [1, 2, 3] -> 3.54 + [0, 1, 2] -> 3.60 + [5, 6] -> 3.64 + [2, 4] -> 4.47 + [4, 6] -> 4.74 + [4, 5, 6] -> 4.74 + [3, 6] -> 5.50 + [2, 6] -> 6.04 + [2, 3, 6] -> 6.04 + [2, 4, 6] -> 6.10 + [0, 4] -> 6.52 + [0, 2, 4] -> 6.52 **Note:** The Delaunay Čech complex can be easily replaced by the :math:`\alpha`-complex, but note that the resulting filtration values will be different. @@ -117,20 +107,18 @@ Then, it is asked to display information about the :math:`\alpha`-complex. .. testcode:: - from gudhi import AlphaComplex - wgt_ac = AlphaComplex(points=[[ 1., -1., -1.], - [-1., 1., -1.], - [-1., -1., 1.], - [ 1., 1., 1.], - [ 2., 2., 2.]], - weights = [4., 4., 4., 4., 1.]) - - stree = wgt_ac.create_simplex_tree() - print('Weighted alpha complex is of dimension ', stree.dimension(), ' - ', - stree.num_simplices(), ' simplices - ', stree.num_vertices(), ' vertices.') - fmt = '%s -> %.2f' + from gudhi import weighted_alpha_complex + stree = weighted_alpha_complex(points=[[ 1., -1., -1.], + [-1., 1., -1.], + [-1., -1., 1.], + [ 1., 1., 1.], + [ 2., 2., 2.]], + weights = [4., 4., 4., 4., 1.]) + + print(f"Weighted alpha is of dimension {stree.dimension()} - {stree.num_simplices()} simplices - ", + f"{stree.num_vertices()} vertices.") for simplex in stree.get_simplices(): - print(fmt % tuple(simplex)) + print('%s -> %.2f' % tuple(simplex)) The output is: diff --git a/src/python/gudhi/delaunay_complex.pyx b/src/python/gudhi/delaunay_complex.pyx index 8989c75f0b..898dc51f13 100644 --- a/src/python/gudhi/delaunay_complex.pyx +++ b/src/python/gudhi/delaunay_complex.pyx @@ -8,7 +8,8 @@ # # Modification(s): # - 2024/03 Vincent Rouvreau: Renamed AlphaComplex as DelaunayComplex. AlphaComplex inherits from it. -# - 2024/10 Vincent Rouvreau: Add square root filtration values interface +# - 2024/10 Vincent Rouvreau: Add square root filtration values interface and new function interfaces (instead of +# classes) # - YYYY/MM Author: Description of the modification from __future__ import print_function @@ -65,12 +66,11 @@ cdef class DelaunayComplex: # Fake constructor that does nothing but documenting the constructor def __init__(self, points: Iterable[Iterable[float]] = [], weights: Optional[Iterable[float]] = None, precision: Literal['fast', 'safe', 'exact'] = 'safe'): - """DelaunayComplex constructor. - - Args: + """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'. + weights (Optional[Iterable[float]]): A list of weights. If set, the number of weights must correspond to + the number of points. + precision (str): Complex precision can be 'fast', 'safe' or 'exact'. Default is 'safe'. :raises ValueError: In case of inconsistency between the number of points and weights. """ @@ -78,7 +78,11 @@ cdef class DelaunayComplex: # The real cython constructor 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'" + 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' @@ -116,9 +120,9 @@ cdef class DelaunayComplex: to `cech` to compute the Delaunay Cech complex. output_squared_values: 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, ...). + 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\'") @@ -140,12 +144,13 @@ cdef class DelaunayComplex: def set_float_relative_precision(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 + precision: When constructing :func:`~gudhi.delaunay_cech_complex`, :func:`~gudhi.alpha_complex`, or + :func:`~gudhi.weighted_alpha_complex` with :code:`precision = 'safe'` (the default), one can + set the float relative precision of filtration values computed. 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 `_ + + :raises ValueError: If precision is not in (0, 1). """ if precision <=0. or precision >= 1.: raise ValueError("Relative precision value must be strictly greater than 0 and strictly lower than 1") @@ -155,8 +160,7 @@ cdef class DelaunayComplex: def get_float_relative_precision() -> float: """ Returns: - The float relative precision of filtration values computation in - :func:`~gudhi.DelaunayComplex.create_simplex_tree` when the DelaunayComplex is constructed with + The float relative precision of filtration values computation when constructing with :code:`precision = 'safe'` (the default). """ return Delaunay_complex_interface.get_float_relative_precision() @@ -176,7 +180,7 @@ cdef class AlphaComplex(DelaunayComplex): .. note:: - When DelaunayComplex is constructed with an infinite value of alpha, the complex is a Delaunay complex. + When Alpha Complex is constructed with an infinite value of alpha, the complex is a Delaunay complex. """ def create_simplex_tree(self, double max_alpha_square: float = float('inf'), default_filtration_value: bool = False, @@ -217,37 +221,132 @@ cdef class AlphaComplex(DelaunayComplex): return self.this_ptr.get_point(vertex) -cdef class DelaunayCechComplex(DelaunayComplex): - """DelaunayCechComplex is a simplicial complex constructed from the finite cells of a Delaunay Triangulation. +def delaunay_cech_complex( + points: Iterable[Iterable[float]] = [], + precision: Literal["fast", "safe", "exact"] = "safe", + max_alpha: float = float("inf"), + output_squared_values: bool = True, +) -> SimplexTree: + """Delaunay Čech complex is a simplicial complex constructed from the finite cells of a Delaunay Triangulation. - The filtration value of each simplex is equal to the squared radius of its minimal enclosing ball (MEB). + The filtration value of each simplex is equal to the radius (squared if `output_squared_values` is set to False) of + its minimal enclosing ball (MEB). - All simplices that have a filtration value strictly greater than a given alpha squared value are not inserted into - the complex. + All simplices that have a filtration value strictly greater than a given alpha value are not inserted into the + complex. .. note:: - When DelaunayCechComplex is constructed with an infinite value of alpha, the complex is a Delaunay complex. + When Delaunay Čech complex is constructed with an infinite value of alpha, the complex is a Delaunay complex. + + 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'. + max_alpha: 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. + output_squared_values: Square root filtration values when True. Default is True, but computation is faster when + set to 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, ...). + """ + cpx = DelaunayComplex(points=points, weights=None, precision=precision) + return cpx.create_simplex_tree( + max_alpha_square=max_alpha * max_alpha, filtration="cech", output_squared_values=output_squared_values + ) + +def delaunay_complex( + points: Iterable[Iterable[float]] = [], +) -> SimplexTree: + """Delaunay Complex is a simplicial complex constructed from the finite cells of a Delaunay Triangulation. This is the method candidate if filtration values are not needed to be computed (will be set to `NaN`). + + Args: + points (Iterable[Iterable[float]]): A list of points in d-Dimension. + + 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, ...). """ - def __init__(self, points: Iterable[Iterable[float]] = [], precision: Literal['fast', 'safe', 'exact'] = 'safe'): - """DelaunayCechComplex constructor. + # Weights and precision has no effect on the triangulation, set precision to "fast" for fatser computation + cpx = DelaunayComplex(points=points, weights=None, precision="fast") + return cpx.create_simplex_tree( + max_alpha_square=0., filtration=None, output_squared_values=False + ) + +def alpha_complex( + points: Iterable[Iterable[float]] = [], + precision: Literal["fast", "safe", "exact"] = "safe", + double max_alpha: float = float("inf"), + bool output_squared_values: bool = True +) -> SimplexTree: + """Alpha complex is a simplicial complex constructed from the finite cells of a Delaunay Triangulation. + + The filtration value of each simplex is computed as the radius of the smallest empty sphere passing through all of + its vertices. + + All simplices that have a filtration value strictly greater than a given alpha value are not inserted into the + complex. - 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) + For more details about the algorithm, please refer to the + `Alpha complex C++ documentation `_ - def create_simplex_tree(self, double max_alpha_square: float = float('inf'), - bool output_squared_values: bool = False) -> 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. - output_squared_values: 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, ...). - """ - return super().create_simplex_tree(max_alpha_square, 'cech', output_squared_values) + .. note:: + + When DelaunayComplex is constructed with an infinite value of alpha, the complex is a Delaunay complex. + + Args: + points (Iterable[Iterable[float]]): A list of points in d-Dimension. + precision (str): Alpha complex precision can be 'fast', 'safe' or 'exact'. Default is 'safe'. + max_alpha: The maximum alpha 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. + output_squared_values: Square root filtration values when True. Default is True, but computation is faster when + set to 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, ...). + """ + cpx = DelaunayComplex(points=points, weights=None, precision=precision) + return cpx.create_simplex_tree( + max_alpha_square=max_alpha * max_alpha, filtration="alpha", output_squared_values=output_squared_values + ) + +def weighted_alpha_complex( + points: Iterable[Iterable[float]] = [], + weights: Iterable[float] = [], + precision: Literal["fast", "safe", "exact"] = "safe", + double max_power_distance: float = float("inf"), +) -> SimplexTree: + """Weighted Alpha complex is a simplicial complex constructed from the finite cells of a Delaunay Triangulation. + + The filtration value of each simplex is computed as the power distance of the smallest power sphere passing through + all of its vertices. + + All simplices that have a filtration value strictly greater than a given alpha square value are not inserted into + the complex. + + For more details about the algorithm, please refer to the + `Alpha complex C++ documentation `_ + + .. note:: + + When DelaunayComplex is constructed with an infinite value of power distance, the complex is a Delaunay complex. + + Args: + points (Iterable[Iterable[float]]): A list of points in d-Dimension. + precision (str): Weighted Alpha complex precision can be 'fast', 'safe' or 'exact'. Default is 'safe'. + 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. + 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, ...). + + :raises ValueError: In case of inconsistency between the number of points and weights. + """ + cpx = DelaunayComplex(points=points, weights=weights, precision=precision) + return cpx.create_simplex_tree( + max_alpha_square=max_power_distance, filtration="alpha", output_squared_values=False + ) diff --git a/src/python/include/Delaunay_complex_factory.h b/src/python/include/Delaunay_complex_factory.h index c945d1637d..1f1df5b475 100644 --- a/src/python/include/Delaunay_complex_factory.h +++ b/src/python/include/Delaunay_complex_factory.h @@ -27,6 +27,7 @@ #include #include // for std::unique_ptr #include +#include // for std::sqrt namespace Gudhi { @@ -85,11 +86,13 @@ 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 - if (output_squared_values) + if (output_squared_values) { Gudhi::cech_complex::assign_MEB_filtration(Kernel(), *simplex_tree, points); - else + simplex_tree->prune_above_filtration(std::sqrt(max_alpha_square)); + } else { Gudhi::cech_complex::assign_MEB_filtration(Kernel(), *simplex_tree, points); - simplex_tree->prune_above_filtration(max_alpha_square); + simplex_tree->prune_above_filtration(max_alpha_square); + } } return result; } else { diff --git a/src/python/test/test_delaunay_complex.py b/src/python/test/test_delaunay_complex.py index 84281fcd11..45a688d437 100644 --- a/src/python/test/test_delaunay_complex.py +++ b/src/python/test/test_delaunay_complex.py @@ -9,7 +9,7 @@ - YYYY/MM Author: Description of the modification """ -from gudhi import AlphaComplex, DelaunayComplex, DelaunayCechComplex +from gudhi import weighted_alpha_complex, alpha_complex, delaunay_complex, delaunay_cech_complex import math import numpy as np import pytest @@ -23,63 +23,97 @@ -def _empty_complex(simplicial_complex, precision): - cplx = simplicial_complex(precision = precision) - assert cplx._is_defined() == True +def test_empty_complex(): + # Specific for delaunay_complex as no precision + stree = delaunay_complex() + assert stree.is_empty() -def _one_2d_point_complex(simplicial_complex, precision): - cplx = simplicial_complex(points=[[0, 0]], precision = precision) - assert cplx._is_defined() == True + stree = delaunay_complex(points=[[0, 0]]) + assert stree.num_vertices() == 1 + assert stree.num_simplices() == 1 -def test_empty_complex(): - for simplicial_complex in [AlphaComplex, DelaunayComplex, DelaunayCechComplex]: - for precision in ['fast', 'safe', 'exact']: - _empty_complex(simplicial_complex, precision) - _one_2d_point_complex(simplicial_complex, precision) + for precision in ['fast', 'safe', 'exact']: + # Specific for delaunay_complex as it requires weights + stree = weighted_alpha_complex(precision = precision) + assert stree.is_empty() + + stree = weighted_alpha_complex(points=[[0, 0]], weights=[0], precision = precision) + assert stree.num_vertices() == 1 + assert stree.num_simplices() == 1 + + for simplicial_complex_helper in [alpha_complex, delaunay_cech_complex]: + stree = simplicial_complex_helper(precision = precision) + assert stree.is_empty() -def _infinite_threshold(simplicial_complex, precision): + stree = simplicial_complex_helper(points=[[0, 0]], precision = precision) + assert stree.num_vertices() == 1 + assert stree.num_simplices() == 1 + + +def _infinite_threshold(simplicial_complex_helper, precision): point_list = [[0, 0], [1, 0], [0, 1], [1, 1]] - cplx = simplicial_complex(points=point_list, precision = precision) - assert cplx._is_defined() == True + simplex_tree = simplicial_complex_helper(points=point_list, precision=precision) - simplex_tree = cplx.create_simplex_tree() assert simplex_tree._is_persistence_defined() == False assert simplex_tree.num_simplices() == 11 assert simplex_tree.num_vertices() == 4 - assert list(simplex_tree.get_filtration()) == [ - ([0], 0.0), - ([1], 0.0), - ([2], 0.0), - ([3], 0.0), - ([0, 1], 0.25), - ([0, 2], 0.25), - ([1, 3], 0.25), - ([2, 3], 0.25), - ([1, 2], 0.5), - ([0, 1, 2], 0.5), - ([1, 2, 3], 0.5), + diag_filt = 1. / math.sqrt(2.) + simplices = [filt[0] for filt in simplex_tree.get_filtration()] + assert simplices == [ + [0], + [1], + [2], + [3], + [0, 1], + [0, 2], + [1, 3], + [2, 3], + [1, 2], + [0, 1, 2], + [1, 2, 3], ] - - assert simplex_tree.get_star([0]) == [ - ([0], 0.0), - ([0, 1], 0.25), - ([0, 1, 2], 0.5), - ([0, 2], 0.25), + filtrations = [filt[1] for filt in simplex_tree.get_filtration()] + np.testing.assert_array_almost_equal(filtrations, [ + 0.0, + 0.0, + 0.0, + 0.0, + 0.5, + 0.5, + 0.5, + 0.5, + diag_filt, + diag_filt, + diag_filt, + ]) + + simplices = [filt[0] for filt in simplex_tree.get_star([0])] + assert simplices == [ + [0], + [0, 1], + [0, 1, 2], + [0, 2], ] - assert simplex_tree.get_cofaces([0], 1) == [([0, 1], 0.25), ([0, 2], 0.25)] + filtrations = [filt[1] for filt in simplex_tree.get_star([0])] + np.testing.assert_array_almost_equal(filtrations, [ + 0.0, + 0.5, + diag_filt, + 0.5, + ]) + + assert simplex_tree.get_cofaces([0], 1) == [([0, 1], 0.5), ([0, 2], 0.5)] def test_infinite_threshold(): - for simplicial_complex in [AlphaComplex, DelaunayCechComplex]: + for simplicial_complex_helper in [alpha_complex, delaunay_cech_complex]: for precision in ['fast', 'safe', 'exact']: - _infinite_threshold(simplicial_complex, precision) + _infinite_threshold(simplicial_complex_helper, precision) -def _finite_threshold(simplicial_complex, precision): +def _finite_threshold(simplicial_complex_helper, precision): point_list = [[0, 0], [1, 0], [0, 1], [1, 1]] - cplx = simplicial_complex(points=point_list, precision = precision) - - simplex_tree = cplx.create_simplex_tree(max_alpha_square=0.25) + simplex_tree = simplicial_complex_helper(points=point_list, max_alpha=0.5, precision=precision) assert simplex_tree.num_simplices() == 8 assert simplex_tree.num_vertices() == 4 @@ -89,20 +123,20 @@ def _finite_threshold(simplicial_complex, precision): ([1], 0.0), ([2], 0.0), ([3], 0.0), - ([0, 1], 0.25), - ([0, 2], 0.25), - ([1, 3], 0.25), - ([2, 3], 0.25), + ([0, 1], 0.5), + ([0, 2], 0.5), + ([1, 3], 0.5), + ([2, 3], 0.5), ] - assert simplex_tree.get_star([0]) == [([0], 0.0), ([0, 1], 0.25), ([0, 2], 0.25)] - assert simplex_tree.get_cofaces([0], 1) == [([0, 1], 0.25), ([0, 2], 0.25)] + assert simplex_tree.get_star([0]) == [([0], 0.0), ([0, 1], 0.5), ([0, 2], 0.5)] + assert simplex_tree.get_cofaces([0], 1) == [([0, 1], 0.5), ([0, 2], 0.5)] def test_filtered_complex(): - for simplicial_complex in [AlphaComplex, DelaunayCechComplex]: + for simplicial_complex_helper in [alpha_complex, delaunay_cech_complex]: for precision in ['fast', 'safe', 'exact']: - _finite_threshold(simplicial_complex, precision) + _finite_threshold(simplicial_complex_helper, precision) -def _safe_persistence_comparison(simplicial_complex, precision): +def _safe_persistence_comparison(simplicial_complex_helper, precision): #generate periodic signal time = np.arange(0, 10, 1) signal = [math.sin(x) for x in time] @@ -113,12 +147,10 @@ def _safe_persistence_comparison(simplicial_complex, precision): 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() + #build simplicial_complex_helper and simplex tree + simplex_tree1 = simplicial_complex_helper(points=embedding1, precision = precision) - cplx2 = simplicial_complex(points=embedding2, precision = precision) - simplex_tree2 = cplx2.create_simplex_tree() + simplex_tree2 = simplicial_complex_helper(points=embedding2, precision = precision) diag1 = simplex_tree1.persistence() diag2 = simplex_tree2.persistence() @@ -129,16 +161,14 @@ def _safe_persistence_comparison(simplicial_complex, precision): def test_safe_persistence_comparison(): - for simplicial_complex in [AlphaComplex, DelaunayCechComplex]: + for simplicial_complex_helper in [alpha_complex, delaunay_cech_complex]: # Won't work for 'fast' version - _safe_persistence_comparison(simplicial_complex, 'safe') - _safe_persistence_comparison(simplicial_complex, 'exact') + _safe_persistence_comparison(simplicial_complex_helper, 'safe') + _safe_persistence_comparison(simplicial_complex_helper, 'exact') -def _delaunay_complex(precision): +def test_delaunay_complex(): point_list = [[0, 0], [1, 0], [0, 1], [1, 1]] - cplx = DelaunayComplex(points=point_list, precision = precision) - - simplex_tree = cplx.create_simplex_tree() + simplex_tree = delaunay_complex(points=point_list) assert simplex_tree.num_simplices() == 11 assert simplex_tree.num_vertices() == 4 @@ -150,56 +180,48 @@ def _delaunay_complex(precision): for filtered_value in simplex_tree.get_cofaces([0], 1): assert math.isnan(filtered_value[1]) -def test_delaunay_complex(): - for precision in ['fast', 'safe', 'exact']: - _delaunay_complex(precision) - -def _3d_points_on_a_plane(simplicial_complex, precision): - cplx = simplicial_complex(points = [[1.0, 1.0 , 0.0], - [7.0, 0.0 , 0.0], - [4.0, 6.0 , 0.0], - [9.0, 6.0 , 0.0], - [0.0, 14.0, 0.0], - [2.0, 19.0, 0.0], - [9.0, 17.0, 0.0]], precision = precision) - - simplex_tree = cplx.create_simplex_tree() +def _3d_points_on_a_plane(simplicial_complex_helper): + simplex_tree = simplicial_complex_helper(points = [[1.0, 1.0 , 0.0], + [7.0, 0.0 , 0.0], + [4.0, 6.0 , 0.0], + [9.0, 6.0 , 0.0], + [0.0, 14.0, 0.0], + [2.0, 19.0, 0.0], + [9.0, 17.0, 0.0]]) + assert simplex_tree.dimension() == 2 assert simplex_tree.num_vertices() == 7 assert simplex_tree.num_simplices() == 25 def test_3d_points_on_a_plane(): - for simplicial_complex in [AlphaComplex, DelaunayComplex, DelaunayCechComplex]: + for simplicial_complex_helper in [alpha_complex, delaunay_complex, delaunay_cech_complex]: for precision in ['fast', 'safe', 'exact']: - _3d_points_on_a_plane(simplicial_complex, precision) - -def _duplicated_2d_points_on_a_plane(simplicial_complex, precision): - cplx = simplicial_complex(points = [[1.0, 1.0 ], - [7.0, 0.0 ], # This point is duplicate - [4.0, 6.0 ], - [9.0, 6.0 ], - [0.0, 14.0], - [2.0, 19.0], - [7.0, 0.0 ], # This point is duplicate - [9.0, 17.0]], precision = precision) - - simplex_tree = cplx.create_simplex_tree() + _3d_points_on_a_plane(simplicial_complex_helper) + +def _duplicated_2d_points_on_a_plane(simplicial_complex_helper): + simplex_tree = simplicial_complex_helper(points = [[1.0, 1.0 ], + [7.0, 0.0 ], # This point is duplicate + [4.0, 6.0 ], + [9.0, 6.0 ], + [0.0, 14.0], + [2.0, 19.0], + [7.0, 0.0 ], # This point is duplicate + [9.0, 17.0]]) + assert simplex_tree.dimension() == 2 assert simplex_tree.num_vertices() == 7 assert simplex_tree.num_simplices() == 25 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) + for simplicial_complex_helper in [alpha_complex, delaunay_complex, delaunay_cech_complex]: + _duplicated_2d_points_on_a_plane(simplicial_complex_helper) def test_output_squared_values(): - for filtration in ['alpha', 'cech', None]: + for simplicial_complex_helper in [alpha_complex, delaunay_cech_complex]: 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', output_squared_values=False) - stree_sqrt = dc.create_simplex_tree(filtration='cech', output_squared_values=True) + pts=[[1, 1], [7, 0], [4, 6], [9, 6], [0, 14], [2, 19], [9, 17]] + stree = simplicial_complex_helper(points=pts, precision=precision, output_squared_values=False) + stree_sqrt = simplicial_complex_helper(points=pts, precision=precision, output_squared_values=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 From e74f9e51ccb1353eaf16f2f901a92a4b4b7937c0 Mon Sep 17 00:00:00 2001 From: Vincent Rouvreau Date: Tue, 26 Nov 2024 12:31:47 +0100 Subject: [PATCH 12/17] Fix some tests and improve test_output_squared_values --- src/python/doc/delaunay_complex_user.rst | 112 +++++++++++------------ src/python/test/test_delaunay_complex.py | 18 ++-- 2 files changed, 67 insertions(+), 63 deletions(-) diff --git a/src/python/doc/delaunay_complex_user.rst b/src/python/doc/delaunay_complex_user.rst index fb2e966c67..f9c4c4a7aa 100644 --- a/src/python/doc/delaunay_complex_user.rst +++ b/src/python/doc/delaunay_complex_user.rst @@ -63,32 +63,32 @@ The output is: .. testoutput:: - Complex is of dimension 2 - 25 simplices - 7 vertices. - [0] -> 0.00 - [1] -> 0.00 - [2] -> 0.00 - [3] -> 0.00 - [4] -> 0.00 - [5] -> 0.00 - [6] -> 0.00 - [2, 3] -> 2.50 - [4, 5] -> 2.69 - [0, 2] -> 2.92 - [0, 1] -> 3.04 - [1, 3] -> 3.16 - [1, 2] -> 3.35 - [1, 2, 3] -> 3.54 - [0, 1, 2] -> 3.60 - [5, 6] -> 3.64 - [2, 4] -> 4.47 - [4, 6] -> 4.74 - [4, 5, 6] -> 4.74 - [3, 6] -> 5.50 - [2, 6] -> 6.04 - [2, 3, 6] -> 6.04 - [2, 4, 6] -> 6.10 - [0, 4] -> 6.52 - [0, 2, 4] -> 6.52 + Complex is of dimension 2 - 25 simplices - 7 vertices. + [0] -> 0.00 + [1] -> 0.00 + [2] -> 0.00 + [3] -> 0.00 + [4] -> 0.00 + [5] -> 0.00 + [6] -> 0.00 + [2, 3] -> 2.50 + [4, 5] -> 2.69 + [0, 2] -> 2.92 + [0, 1] -> 3.04 + [1, 3] -> 3.16 + [1, 2] -> 3.35 + [1, 2, 3] -> 3.54 + [0, 1, 2] -> 3.60 + [5, 6] -> 3.64 + [2, 4] -> 4.47 + [4, 6] -> 4.74 + [4, 5, 6] -> 4.74 + [3, 6] -> 5.50 + [2, 6] -> 6.04 + [2, 3, 6] -> 6.04 + [2, 4, 6] -> 6.10 + [0, 4] -> 6.52 + [0, 2, 4] -> 6.52 **Note:** The Delaunay Čech complex can be easily replaced by the :math:`\alpha`-complex, but note that the resulting filtration values will be different. @@ -124,33 +124,33 @@ The output is: .. testoutput:: - Weighted alpha complex is of dimension 3 - 29 simplices - 5 vertices. - [0, 1, 2, 3] -> -1.00 - [0, 1, 2] -> -1.33 - [0, 1, 3, 4] -> 95.00 - [0, 1, 3] -> -1.33 - [0, 1, 4] -> 95.00 - [0, 1] -> -2.00 - [0, 2, 3, 4] -> 95.00 - [0, 2, 3] -> -1.33 - [0, 2, 4] -> 95.00 - [0, 2] -> -2.00 - [0, 3, 4] -> 23.00 - [0, 3] -> -2.00 - [0, 4] -> 23.00 - [0] -> -4.00 - [1, 2, 3, 4] -> 95.00 - [1, 2, 3] -> -1.33 - [1, 2, 4] -> 95.00 - [1, 2] -> -2.00 - [1, 3, 4] -> 23.00 - [1, 3] -> -2.00 - [1, 4] -> 23.00 - [1] -> -4.00 - [2, 3, 4] -> 23.00 - [2, 3] -> -2.00 - [2, 4] -> 23.00 - [2] -> -4.00 - [3, 4] -> -1.00 - [3] -> -4.00 - [4] -> -1.00 + Weighted alpha is of dimension 3 - 29 simplices - 5 vertices. + [0, 1, 2, 3] -> -1.00 + [0, 1, 2] -> -1.33 + [0, 1, 3, 4] -> 95.00 + [0, 1, 3] -> -1.33 + [0, 1, 4] -> 95.00 + [0, 1] -> -2.00 + [0, 2, 3, 4] -> 95.00 + [0, 2, 3] -> -1.33 + [0, 2, 4] -> 95.00 + [0, 2] -> -2.00 + [0, 3, 4] -> 23.00 + [0, 3] -> -2.00 + [0, 4] -> 23.00 + [0] -> -4.00 + [1, 2, 3, 4] -> 95.00 + [1, 2, 3] -> -1.33 + [1, 2, 4] -> 95.00 + [1, 2] -> -2.00 + [1, 3, 4] -> 23.00 + [1, 3] -> -2.00 + [1, 4] -> 23.00 + [1] -> -4.00 + [2, 3, 4] -> 23.00 + [2, 3] -> -2.00 + [2, 4] -> 23.00 + [2] -> -4.00 + [3, 4] -> -1.00 + [3] -> -4.00 + [4] -> -1.00 diff --git a/src/python/test/test_delaunay_complex.py b/src/python/test/test_delaunay_complex.py index 45a688d437..066d14d7a4 100644 --- a/src/python/test/test_delaunay_complex.py +++ b/src/python/test/test_delaunay_complex.py @@ -219,10 +219,14 @@ def test_duplicated_2d_points_on_a_plane(): def test_output_squared_values(): for simplicial_complex_helper in [alpha_complex, delaunay_cech_complex]: for precision in ['fast', 'safe', 'exact']: - pts=[[1, 1], [7, 0], [4, 6], [9, 6], [0, 14], [2, 19], [9, 17]] - stree = simplicial_complex_helper(points=pts, precision=precision, output_squared_values=False) - stree_sqrt = simplicial_complex_helper(points=pts, precision=precision, output_squared_values=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))) + for max_alpha in [float('inf'), math.sqrt(20.)]: + pts=[[1, 1], [7, 0], [4, 6], [9, 6], [0, 14], [2, 19], [9, 17]] + stree = simplicial_complex_helper(points=pts, precision=precision, + output_squared_values=False, max_alpha=max_alpha) + stree_sqrt = simplicial_complex_helper(points=pts, precision=precision, + output_squared_values=True, max_alpha=max_alpha) + assert stree.num_simplices() == stree_sqrt.num_simplices() + 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))) From b2c979400f54bca2b3026e124800d3ca9ed84e6e Mon Sep 17 00:00:00 2001 From: Vincent Rouvreau Date: Tue, 28 Jan 2025 09:23:18 +0100 Subject: [PATCH 13/17] code review: Remove duplicated type mypy+docstring, mypy is sufficient enough --- src/python/gudhi/delaunay_complex.pyx | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) diff --git a/src/python/gudhi/delaunay_complex.pyx b/src/python/gudhi/delaunay_complex.pyx index 898dc51f13..c797461bfb 100644 --- a/src/python/gudhi/delaunay_complex.pyx +++ b/src/python/gudhi/delaunay_complex.pyx @@ -240,8 +240,8 @@ def delaunay_cech_complex( When Delaunay Čech complex is constructed with an infinite value of alpha, the complex is a Delaunay complex. 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'. + points: A list of points in d-Dimension. + precision: Delaunay Čech complex precision can be 'fast', 'safe' or 'exact'. Default is 'safe'. max_alpha: 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. output_squared_values: Square root filtration values when True. Default is True, but computation is faster when @@ -262,7 +262,7 @@ def delaunay_complex( """Delaunay Complex is a simplicial complex constructed from the finite cells of a Delaunay Triangulation. This is the method candidate if filtration values are not needed to be computed (will be set to `NaN`). Args: - points (Iterable[Iterable[float]]): A list of points in d-Dimension. + points: A list of points in d-Dimension. Returns: SimplexTree: A simplex tree created from the Delaunay Triangulation. The vertex `k` corresponds to the k-th @@ -297,8 +297,8 @@ def alpha_complex( When DelaunayComplex is constructed with an infinite value of alpha, the complex is a Delaunay complex. Args: - points (Iterable[Iterable[float]]): A list of points in d-Dimension. - precision (str): Alpha complex precision can be 'fast', 'safe' or 'exact'. Default is 'safe'. + points: A list of points in d-Dimension. + precision: Alpha complex precision can be 'fast', 'safe' or 'exact'. Default is 'safe'. max_alpha: The maximum alpha 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. output_squared_values: Square root filtration values when True. Default is True, but computation is faster when @@ -335,9 +335,10 @@ def weighted_alpha_complex( When DelaunayComplex is constructed with an infinite value of power distance, the complex is a Delaunay complex. Args: - points (Iterable[Iterable[float]]): A list of points in d-Dimension. - precision (str): Weighted Alpha complex precision can be 'fast', 'safe' or 'exact'. Default is 'safe'. - max_alpha_square: The maximum alpha square threshold the simplices shall not exceed. Default is set to + points: A list of points in d-Dimension. + weights: A list of weights. If set, the number of weights must correspond to the number of points. + precision: Weighted Alpha complex precision can be 'fast', 'safe' or 'exact'. Default is 'safe'. + max_power_distance: 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. Returns: SimplexTree: A simplex tree created from the Delaunay Triangulation. The vertex `k` corresponds to the k-th From 4ab7e4ece648fff8de9dc70fb2cae62703f7dad1 Mon Sep 17 00:00:00 2001 From: Vincent Rouvreau Date: Tue, 28 Jan 2025 16:54:55 +0100 Subject: [PATCH 14/17] code review: confusion true/false when rename + doc for Alpha complex --- src/Alpha_complex/include/gudhi/Alpha_complex.h | 17 ++++++++++------- .../test/Delaunay_complex_unit_test.h | 3 ++- .../test/Weighted_alpha_complex_unit_test.cpp | 3 ++- .../utilities/alpha_complex_persistence.cpp | 16 ++++++++++------ src/Alpha_complex/utilities/alphacomplex.md | 6 +++--- src/Cech_complex/utilities/cechcomplex.md | 2 +- .../utilities/delaunay_cech_persistence.cpp | 4 ++-- 7 files changed, 30 insertions(+), 21 deletions(-) diff --git a/src/Alpha_complex/include/gudhi/Alpha_complex.h b/src/Alpha_complex/include/gudhi/Alpha_complex.h index e9b9ab8c43..e8ec599034 100644 --- a/src/Alpha_complex/include/gudhi/Alpha_complex.h +++ b/src/Alpha_complex/include/gudhi/Alpha_complex.h @@ -355,12 +355,12 @@ class Alpha_complex { * It also computes the filtration values accordingly to the \ref createcomplexalgorithm if default_filtration_value * is not set. * - * \tparam output_squared_values 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 output_squared_values is true. + * \tparam output_squared_values If `true` (default value), it assigns to each simplex a filtration value equal to + * \f$ \alpha^2 \f$, or to \f$ \alpha \f$ when `output_squared_values` is `false`. cf. \ref createcomplexalgorithm * \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 + * @param[in] max_alpha_square maximum for \f$ \alpha^2 \f$ 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 * `true`. * @param[in] exact Exact filtration values computation. Not exact if `Kernel` is not bool create_complex(SimplicialComplexForAlpha& complex, @@ -464,7 +467,7 @@ class Alpha_complex { if(exact) CGAL::exact(sqrad); #endif alpha_complex_filtration = cgal_converter(sqrad); - if constexpr (output_squared_values) { + if constexpr (!output_squared_values) { alpha_complex_filtration = sqrt(alpha_complex_filtration); } } @@ -489,7 +492,7 @@ class Alpha_complex { // 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 - if constexpr (output_squared_values) { + if constexpr (!output_squared_values) { complex.prune_above_filtration(sqrt(max_alpha_square)); } else { complex.prune_above_filtration(max_alpha_square); diff --git a/src/Alpha_complex/test/Delaunay_complex_unit_test.h b/src/Alpha_complex/test/Delaunay_complex_unit_test.h index aa16423b1a..e58cc1e4eb 100644 --- a/src/Alpha_complex/test/Delaunay_complex_unit_test.h +++ b/src/Alpha_complex/test/Delaunay_complex_unit_test.h @@ -49,7 +49,8 @@ void compare_delaunay_complex_simplices() { 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)); + // set output_squared_values to false + BOOST_CHECK(alpha_complex.template create_complex(stree_from_alpha_sqrt)); std::clog << "Check simplices from alpha complex filtration values when output_squared_values is true\n"; // Check all the simplices from alpha complex are in the Delaunay complex 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 8a353228fc..120906b564 100644 --- a/src/Alpha_complex/test/Weighted_alpha_complex_unit_test.cpp +++ b/src/Alpha_complex/test/Weighted_alpha_complex_unit_test.cpp @@ -156,7 +156,8 @@ BOOST_AUTO_TEST_CASE(Is_weighted_alpha_complex_nan) { } std::clog << "Weighted alpha complex with output_squared_values\n"; Gudhi::Simplex_tree<> stree_sqrt; - if (alpha_complex_from_weighted_points.create_complex(stree_sqrt)) { + // set output_squared_values to false means that negative squared values will be NaN + 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)) { diff --git a/src/Alpha_complex/utilities/alpha_complex_persistence.cpp b/src/Alpha_complex/utilities/alpha_complex_persistence.cpp index b8cc6fc19a..1c0cbb13fa 100644 --- a/src/Alpha_complex/utilities/alpha_complex_persistence.cpp +++ b/src/Alpha_complex/utilities/alpha_complex_persistence.cpp @@ -80,17 +80,21 @@ Simplex_tree create_simplex_tree(const std::string &off_file_points, const std:: Gudhi::alpha_complex::Alpha_complex alpha_complex_from_file(points, weights); if (output_squared_values) - complex_creation = alpha_complex_from_file.template create_complex(stree, alpha_square_max_value, exact_version); + 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); + complex_creation = alpha_complex_from_file.template 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 (output_squared_values) - complex_creation = alpha_complex_from_file.template create_complex(stree, alpha_square_max_value, exact_version); + 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); + complex_creation = alpha_complex_from_file.template create_complex(stree, alpha_square_max_value, + exact_version); } if (!complex_creation) { std::cerr << "Alpha complex simplicial complex creation failed." << std::endl; @@ -205,9 +209,9 @@ void program_options(int argc, char *argv[], std::string &off_file_points, bool 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 << "'--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 << "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"; diff --git a/src/Alpha_complex/utilities/alphacomplex.md b/src/Alpha_complex/utilities/alphacomplex.md index 5272b355e3..af358c2509 100644 --- a/src/Alpha_complex/utilities/alphacomplex.md +++ b/src/Alpha_complex/utilities/alphacomplex.md @@ -22,10 +22,10 @@ Different versions of Alpha complex computation are available: 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. +'--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 +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: @@ -68,7 +68,7 @@ points (one value per line). Default version is not weighted. **Example** ``` - alpha_complex_persistence -r 32 -p 2 -m 0.45 ../../data/points/tore3D_300.off + alpha_complex_persistence -p 2 -m 0.45 ../../data/points/tore3D_300.off ``` N.B.: diff --git a/src/Cech_complex/utilities/cechcomplex.md b/src/Cech_complex/utilities/cechcomplex.md index 290650d9df..005c7baa10 100644 --- a/src/Cech_complex/utilities/cechcomplex.md +++ b/src/Cech_complex/utilities/cechcomplex.md @@ -72,7 +72,7 @@ Different versions of Delaunay-Čech complex computation are available: 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' +If you are interested in radius of the MEB as filtration values, pass the '--square-root-filtrations' (or '-s') option. diff --git a/src/Cech_complex/utilities/delaunay_cech_persistence.cpp b/src/Cech_complex/utilities/delaunay_cech_persistence.cpp index fa20854dea..c377451187 100644 --- a/src/Cech_complex/utilities/delaunay_cech_persistence.cpp +++ b/src/Cech_complex/utilities/delaunay_cech_persistence.cpp @@ -163,8 +163,8 @@ void program_options(int argc, char* argv[], std::string& off_file_points, bool& 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 << "If you are interested in radius of the MEB as filtration values, pass the "; + std< Date: Tue, 28 Jan 2025 17:55:14 +0100 Subject: [PATCH 15/17] code review: confusion true/false when rename + doc for Cech complex --- src/Cech_complex/include/gudhi/MEB_filtration.h | 14 +++++++------- src/Cech_complex/test/test_cech_complex.cpp | 12 +++++++----- .../utilities/delaunay_cech_persistence.cpp | 2 +- 3 files changed, 15 insertions(+), 13 deletions(-) diff --git a/src/Cech_complex/include/gudhi/MEB_filtration.h b/src/Cech_complex/include/gudhi/MEB_filtration.h index 6f61f68664..013b1362ab 100644 --- a/src/Cech_complex/include/gudhi/MEB_filtration.h +++ b/src/Cech_complex/include/gudhi/MEB_filtration.h @@ -30,8 +30,8 @@ namespace Gudhi::cech_complex { * Applied on a Čech complex, it recomputes the same values (squared or not in function of `output_squared_values`). * Applied on a Delaunay triangulation, it computes the Delaunay-Čech filtration. * - * \tparam output_squared_values 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 output_squared_values is true. + * \tparam output_squared_values If `true` (default value), it assigns to each simplex a filtration value equal to + * the squared radius of the MEB, or to the radius when `output_squared_values` is `false`. * \tparam Kernel CGAL kernel: either Epick_d or Epeck_d. * \tparam PointRange Random access range of `Kernel::Point_d`. * @@ -42,7 +42,7 @@ namespace Gudhi::cech_complex { * 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; @@ -82,7 +82,7 @@ void assign_MEB_filtration(Kernel&&k, SimplicialComplexForMEB& complex, PointRan if (exact) CGAL::exact(r); complex.assign_key(sh, cache_.size()); Filtration_value filt{max(cvt(r), Filtration_value(0))}; - if constexpr (output_squared_values) + if constexpr (!output_squared_values) filt = sqrt(filt); complex.assign_filtration(sh, filt); cache_.emplace_back(std::move(m), std::move(r)); @@ -121,9 +121,9 @@ 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. - Filtration_value filt{max(cvt(r), Filtration_value(0))}; - if constexpr (output_squared_values) - filt = sqrt(filt); + Filtration_value filt{cvt(r)}; + if constexpr (!output_squared_values) + filt = sqrt(max(filt, Filtration_value(0))); // maxf = filt except for rounding errors maxf = max(maxf, filt); complex.assign_key(sh, cache_.size()); diff --git a/src/Cech_complex/test/test_cech_complex.cpp b/src/Cech_complex/test/test_cech_complex.cpp index 061355d6ce..80a50b8b09 100644 --- a/src/Cech_complex/test/test_cech_complex.cpp +++ b/src/Cech_complex/test/test_cech_complex.cpp @@ -182,9 +182,11 @@ BOOST_AUTO_TEST_CASE(Cech_complex_for_documentation) { } st3.reset_filtration(-1); // unnecessary, but ensures we don't cheat + // Means output_squared_values=true 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))); + // sqrt all filtration values + st3.for_each_simplex([&](auto sh, int){ + 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 << " ( "; @@ -193,10 +195,10 @@ BOOST_AUTO_TEST_CASE(Cech_complex_for_documentation) { } BOOST_CHECK(st3 == st3_save); // Should only be an approximate test - // Same test but with output_squared_values=true + // Same test but with output_squared_values=false 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 output_squared_values=true\n"; + Gudhi::cech_complex::assign_MEB_filtration(Kernel(), st3, points); + std::clog << "Simplex tree from assign_MEB_filtration with output_squared_values=false\n"; for (auto f_simplex : st3.filtration_simplex_range()) { std::clog << " ( "; for (auto vertex : st3.simplex_vertex_range(f_simplex)) std::clog << vertex << " "; diff --git a/src/Cech_complex/utilities/delaunay_cech_persistence.cpp b/src/Cech_complex/utilities/delaunay_cech_persistence.cpp index c377451187..64a51a3b3c 100644 --- a/src/Cech_complex/utilities/delaunay_cech_persistence.cpp +++ b/src/Cech_complex/utilities/delaunay_cech_persistence.cpp @@ -164,7 +164,7 @@ void program_options(int argc, char* argv[], std::string& off_file_points, bool& 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 "; - std< Date: Wed, 29 Jan 2025 15:12:01 +0100 Subject: [PATCH 16/17] code review: confusion true/false when rename + doc for python module --- src/python/gudhi/delaunay_complex.pyx | 47 ++++++++++--------- src/python/include/Delaunay_complex_factory.h | 6 +-- src/python/test/test_delaunay_complex.py | 4 +- 3 files changed, 30 insertions(+), 27 deletions(-) diff --git a/src/python/gudhi/delaunay_complex.pyx b/src/python/gudhi/delaunay_complex.pyx index c797461bfb..52dfc80d30 100644 --- a/src/python/gudhi/delaunay_complex.pyx +++ b/src/python/gudhi/delaunay_complex.pyx @@ -70,7 +70,7 @@ cdef class DelaunayComplex: 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): Complex precision can be 'fast', 'safe' or 'exact'. Default is 'safe'. + precision (str): Complex precision can be `'fast'`, `'safe'` or `'exact'`. Default is `'safe'`. :raises ValueError: In case of inconsistency between the number of points and weights. """ @@ -104,7 +104,7 @@ cdef class DelaunayComplex: del self.this_ptr def _is_defined(self): - """Returns true if DelaunayComplex pointer is not NULL. + """Returns `True` if DelaunayComplex pointer is not NULL. """ return self.this_ptr != NULL @@ -118,7 +118,7 @@ cdef class DelaunayComplex: 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. - output_squared_values: Square root filtration values when True. Default is False. + output_squared_values: 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 @@ -184,7 +184,7 @@ cdef class AlphaComplex(DelaunayComplex): """ def create_simplex_tree(self, double max_alpha_square: float = float('inf'), default_filtration_value: bool = False, - bool output_squared_values: bool = False) -> SimplexTree: + bool output_squared_values: bool = True) -> SimplexTree: """ Args: max_alpha_square: The maximum alpha square threshold the simplices shall not exceed. Default is set to @@ -192,7 +192,8 @@ cdef class AlphaComplex(DelaunayComplex): 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. - output_squared_values: Square root filtration values when True. Default is False. + output_squared_values: Square root filtration values when `True`. Default is `True` to keep backward + compatibility. 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 @@ -201,8 +202,8 @@ cdef class AlphaComplex(DelaunayComplex): filtration = 'alpha' if default_filtration_value: filtration = None - warnings.warn('''Since Gudhi 3.10, creating an AlphaComplex with default_filtration_value=True is deprecated. - Please consider constructing a DelaunayComplex instead. + 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, output_squared_values) @@ -225,12 +226,12 @@ def delaunay_cech_complex( points: Iterable[Iterable[float]] = [], precision: Literal["fast", "safe", "exact"] = "safe", max_alpha: float = float("inf"), - output_squared_values: bool = True, + output_squared_values: bool = False, ) -> SimplexTree: """Delaunay Čech complex is a simplicial complex constructed from the finite cells of a Delaunay Triangulation. - The filtration value of each simplex is equal to the radius (squared if `output_squared_values` is set to False) of - its minimal enclosing ball (MEB). + The filtration value of each simplex is equal to the radius (squared if `output_squared_values` is set to `True`) + of its minimal enclosing ball (MEB). All simplices that have a filtration value strictly greater than a given alpha value are not inserted into the complex. @@ -241,15 +242,15 @@ def delaunay_cech_complex( Args: points: A list of points in d-Dimension. - precision: Delaunay Čech complex precision can be 'fast', 'safe' or 'exact'. Default is 'safe'. + precision: Delaunay Čech complex precision can be `'fast'`, `'safe'` or `'exact'`. Default is `'safe'`. max_alpha: 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. - output_squared_values: Square root filtration values when True. Default is True, but computation is faster when - set to False. + output_squared_values: Square root filtration values when `False`. Default is `False`, but computation is + faster when set to `True`. 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, ...). + (duplicate points, ...). """ cpx = DelaunayComplex(points=points, weights=None, precision=precision) return cpx.create_simplex_tree( @@ -279,12 +280,12 @@ def alpha_complex( points: Iterable[Iterable[float]] = [], precision: Literal["fast", "safe", "exact"] = "safe", double max_alpha: float = float("inf"), - bool output_squared_values: bool = True + bool output_squared_values: bool = False ) -> SimplexTree: """Alpha complex is a simplicial complex constructed from the finite cells of a Delaunay Triangulation. - The filtration value of each simplex is computed as the radius of the smallest empty sphere passing through all of - its vertices. + The filtration value of each simplex is computed as the radius (squared if `output_squared_values` is set to + `True`) of the smallest empty sphere passing through all of its vertices. All simplices that have a filtration value strictly greater than a given alpha value are not inserted into the complex. @@ -298,11 +299,11 @@ def alpha_complex( Args: points: A list of points in d-Dimension. - precision: Alpha complex precision can be 'fast', 'safe' or 'exact'. Default is 'safe'. + precision: Alpha complex precision can be `'fast'`, `'safe'` or `'exact'`. Default is `'safe'`. max_alpha: The maximum alpha 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. - output_squared_values: Square root filtration values when True. Default is True, but computation is faster when - set to False. + output_squared_values: Square root filtration values when `False`. Default is `False`, but computation is + faster when set to `True`. 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 @@ -337,7 +338,7 @@ def weighted_alpha_complex( Args: points: A list of points in d-Dimension. weights: A list of weights. If set, the number of weights must correspond to the number of points. - precision: Weighted Alpha complex precision can be 'fast', 'safe' or 'exact'. Default is 'safe'. + precision: Weighted Alpha complex precision can be `'fast'`, `'safe'` or `'exact'`. Default is `'safe'`. max_power_distance: 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. Returns: @@ -348,6 +349,8 @@ def weighted_alpha_complex( :raises ValueError: In case of inconsistency between the number of points and weights. """ cpx = DelaunayComplex(points=points, weights=weights, precision=precision) + # There is no added value to set output_squared_values to False, as filtration values can be negative in a weighted + # case. If output_squared_values is set to False, filtration values would be NaN. return cpx.create_simplex_tree( - max_alpha_square=max_power_distance, filtration="alpha", output_squared_values=False + max_alpha_square=max_power_distance, filtration="alpha", output_squared_values=True ) diff --git a/src/python/include/Delaunay_complex_factory.h b/src/python/include/Delaunay_complex_factory.h index 1f1df5b475..27bdc9e590 100644 --- a/src/python/include/Delaunay_complex_factory.h +++ b/src/python/include/Delaunay_complex_factory.h @@ -86,11 +86,11 @@ 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 - if (output_squared_values) { - Gudhi::cech_complex::assign_MEB_filtration(Kernel(), *simplex_tree, points); + if (!output_squared_values) { + Gudhi::cech_complex::assign_MEB_filtration(Kernel(), *simplex_tree, points); simplex_tree->prune_above_filtration(std::sqrt(max_alpha_square)); } else { - Gudhi::cech_complex::assign_MEB_filtration(Kernel(), *simplex_tree, points); + Gudhi::cech_complex::assign_MEB_filtration(Kernel(), *simplex_tree, points); simplex_tree->prune_above_filtration(max_alpha_square); } } diff --git a/src/python/test/test_delaunay_complex.py b/src/python/test/test_delaunay_complex.py index 066d14d7a4..a7e33e159f 100644 --- a/src/python/test/test_delaunay_complex.py +++ b/src/python/test/test_delaunay_complex.py @@ -222,9 +222,9 @@ def test_output_squared_values(): for max_alpha in [float('inf'), math.sqrt(20.)]: pts=[[1, 1], [7, 0], [4, 6], [9, 6], [0, 14], [2, 19], [9, 17]] stree = simplicial_complex_helper(points=pts, precision=precision, - output_squared_values=False, max_alpha=max_alpha) + output_squared_values=True, max_alpha=max_alpha) stree_sqrt = simplicial_complex_helper(points=pts, precision=precision, - output_squared_values=True, max_alpha=max_alpha) + output_squared_values=False, max_alpha=max_alpha) assert stree.num_simplices() == stree_sqrt.num_simplices() for simplex, filt in stree_sqrt.get_filtration(): # np.testing.assert_almost_equal(float('nan'), float('nan')) is ok From db99ef86c293f4057fc606f4a7773609a2b220e1 Mon Sep 17 00:00:00 2001 From: Vincent Rouvreau Date: Thu, 30 Jan 2025 14:20:41 +0100 Subject: [PATCH 17/17] clean: remove unused import, do not document __new__ --- src/python/doc/delaunay_complex_ref.rst | 1 + src/python/gudhi/delaunay_complex.pyx | 1 - 2 files changed, 1 insertion(+), 1 deletion(-) diff --git a/src/python/doc/delaunay_complex_ref.rst b/src/python/doc/delaunay_complex_ref.rst index fb18ae047a..ead5862045 100644 --- a/src/python/doc/delaunay_complex_ref.rst +++ b/src/python/doc/delaunay_complex_ref.rst @@ -23,6 +23,7 @@ Alpha complex reference manual .. autofunction:: gudhi.weighted_alpha_complex .. autoclass:: gudhi.AlphaComplex + :exclude-members: __new__ .. deprecated:: 3.11.0 Use :py:func:`gudhi.alpha_complex`, :py:func:`gudhi.weighted_alpha_complex` or :py:func:`gudhi.delaunay_complex` diff --git a/src/python/gudhi/delaunay_complex.pyx b/src/python/gudhi/delaunay_complex.pyx index 52dfc80d30..bd4a220209 100644 --- a/src/python/gudhi/delaunay_complex.pyx +++ b/src/python/gudhi/delaunay_complex.pyx @@ -25,7 +25,6 @@ from collections.abc import Iterable from gudhi.simplex_tree cimport * from gudhi.simplex_tree import SimplexTree -from gudhi import read_points_from_off_file __author__ = "Vincent Rouvreau" __copyright__ = "Copyright (C) 2016 Inria"