diff --git a/src/Cech_complex/concept/SimplicialComplexForMEB.h b/src/Cech_complex/concept/SimplicialComplexForMEB.h new file mode 100644 index 0000000000..19a0a90818 --- /dev/null +++ b/src/Cech_complex/concept/SimplicialComplexForMEB.h @@ -0,0 +1,58 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): Marc Glisse + * + * Copyright (C) 2018 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#ifndef CONCEPT_CECH_COMPLEX_SIMPLICIAL_COMPLEX_FOR_MEB_H_ +#define CONCEPT_CECH_COMPLEX_SIMPLICIAL_COMPLEX_FOR_MEB_H_ + +namespace Gudhi { + +namespace cech_complex { + +/** The concept SimplicialComplexForMEB describes the requirements for a type to implement a simplicial + * complex that can be filled by `assign_MEB_filtration()`. It is typically satisfied by `Simplex_tree` if + * `SimplexTreeOptions::store_key` is true. + */ +struct SimplicialComplexForMEB { + /** \brief Handle for a simplex. */ + typedef unspecified Simplex_handle; + /** \brief Handle for a vertex. Must be a non-negative integer, + * it is also used as an index into the input list of points. */ + typedef unspecified Vertex_handle; + /** \brief Type of filtration values. */ + typedef unspecified Filtration_value; + /** \brief Integer type large enough to index all simplices. */ + typedef unspecified Simplex_key; + + /** \brief Returns the filtration value to the 'simplex'. */ + Filtration_value filtration(Simplex_handle simplex); + /** \brief Assigns this 'filtration' value to the 'simplex'. */ + int assign_filtration(Simplex_handle simplex, Filtration_value filtration); + + /** \brief Returns the key assigned to the 'simplex' with `assign_key()`. */ + Simplex_key key(Simplex_handle simplex); + /** \brief Assigns this 'key' to the 'simplex'. */ + void assign_key(Simplex_handle simplex, Simplex_key key); + + /** \brief Returns a range over vertices (as Vertex_handle) of a given simplex. */ + Simplex_vertex_range simplex_vertex_range(Simplex_handle simplex); + + /** \brief Returns a range of the pairs (simplex, opposite vertex) of the boundary of the 'simplex'. */ + Boundary_opposite_vertex_simplex_range boundary_opposite_vertex_simplex_range(Simplex_handle simplex); + + /** \brief Calls `callback(simplex, dim)` for every simplex of the complex, + * with the guarantee that faces are visited before cofaces. */ + void for_each_simplex(auto callback); +}; + +} // namespace cech_complex + +} // namespace Gudhi + +#endif // CONCEPT_CECH_COMPLEX_SIMPLICIAL_COMPLEX_FOR_MEB_H_ diff --git a/src/Cech_complex/doc/Intro_cech_complex.h b/src/Cech_complex/doc/Intro_cech_complex.h index 73093c07d7..c8878d6517 100644 --- a/src/Cech_complex/doc/Intro_cech_complex.h +++ b/src/Cech_complex/doc/Intro_cech_complex.h @@ -17,7 +17,7 @@ namespace cech_complex { /** \defgroup cech_complex Čech complex * - * \author Vincent Rouvreau, Hind montassif + * \author Vincent Rouvreau, Hind montassif, Marc Glisse * * @{ * @@ -62,6 +62,11 @@ namespace cech_complex { * This radius computation is the reason why the Cech_complex is taking much more time to be computed than the * \ref rips_complex but it offers more topological guarantees. * + * If you already have a simplicial complex, it is possible to assign to each simplex a filtration value corresponding + * to the squared radius of its minimal enclosing ball using `assign_MEB_filtration()`. This can provide an alternate + * way of computing a Čech filtration, or it can be used on a Delaunay triangulation to compute a Delaunay-Čech + * filtration. + * * \subsection cechpointscloudexample Example from a point cloud * * This example builds the proximity graph from the given points, and maximal radius values. diff --git a/src/Cech_complex/include/gudhi/Cech_complex_blocker.h b/src/Cech_complex/include/gudhi/Cech_complex_blocker.h index e78e37b768..475def561b 100644 --- a/src/Cech_complex/include/gudhi/Cech_complex_blocker.h +++ b/src/Cech_complex/include/gudhi/Cech_complex_blocker.h @@ -48,7 +48,7 @@ class Cech_blocker { // Numeric type of coordinates in the kernel using FT = typename Kernel::FT; // Sphere is a pair of point and squared radius. - using Sphere = typename std::pair; + using Sphere = std::pair; private: diff --git a/src/Cech_complex/include/gudhi/MEB_filtration.h b/src/Cech_complex/include/gudhi/MEB_filtration.h new file mode 100644 index 0000000000..86d14df86d --- /dev/null +++ b/src/Cech_complex/include/gudhi/MEB_filtration.h @@ -0,0 +1,112 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): Marc Glisse + * + * Copyright (C) 2023 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#ifndef MEB_FILTRATION_H_ +#define MEB_FILTRATION_H_ + +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). + * + * Applied on a Čech complex, it recomputes the same values (squared). Applied on a Delaunay triangulation, it computes the Delaunay-Čech filtration. + * + * \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. + */ + +template +void assign_MEB_filtration(Kernel&&k, SimplicialComplexForMEB& complex, PointRange points, bool exact = false) { + using Point_d = typename Kernel::Point_d; + using FT = typename Kernel::FT; + using Sphere = std::pair; + + using Vertex_handle = typename SimplicialComplexForMEB::Vertex_handle; + using Simplex_handle = typename SimplicialComplexForMEB::Simplex_handle; + using Filtration_value = typename SimplicialComplexForMEB::Filtration_value; + + std::vector cache_; + std::vector pts; + CGAL::NT_converter cvt; + + auto fun = [&](Simplex_handle sh, int dim){ + if (dim == 0) complex.assign_filtration(sh, 0); + else if (dim == 1) { + // For a Simplex_tree, this would be a bit faster, but that's probably negligible + // Vertex_handle u = sh->first; Vertex_handle v = self_siblings(sh)->parent(); + auto verts = complex.simplex_vertex_range(sh); + auto vert_it = verts.begin(); + Vertex_handle u = *vert_it; + Vertex_handle v = *++vert_it; + auto&& pu = points[u]; + Point_d m = k.midpoint_d_object()(pu, points[v]); + FT r = k.squared_distance_d_object()(m, pu); + if (exact) CGAL::exact(r); + complex.assign_key(sh, cache_.size()); + complex.assign_filtration(sh, std::max(cvt(r), Filtration_value(0))); + cache_.emplace_back(std::move(m), std::move(r)); + } else { + Filtration_value maxf = 0; // max filtration of the faces + bool found = false; + using std::max; + for (auto face_opposite_vertex : complex.boundary_opposite_vertex_simplex_range(sh)) { + maxf = max(maxf, complex.filtration(face_opposite_vertex.first)); + if (!found) { + auto key = complex.key(face_opposite_vertex.first); + Sphere const& sph = cache_[key]; + if (k.squared_distance_d_object()(sph.first, points[face_opposite_vertex.second]) > sph.second) continue; + found = true; + complex.assign_key(sh, key); + // With exact computations, we could stop here + // complex.assign_filtration(sh, complex.filtration(face_opposite_vertex.first)); return; + // but because of possible rounding errors, we continue with the equivalent of make_filtration_non_decreasing + } + } + if (!found) { + // None of the faces are good enough, MEB must be the circumsphere. + pts.clear(); + for (auto vertex : complex.simplex_vertex_range(sh)) + pts.push_back(points[vertex]); + Point_d c = k.construct_circumcenter_d_object()(pts.begin(), pts.end()); + FT r = k.squared_distance_d_object()(c, pts.front()); + if (exact) CGAL::exact(r); + maxf = max(maxf, cvt(r)); // maxf = cvt(r) except for rounding errors + 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)); + } + complex.assign_filtration(sh, maxf); + } + }; + complex.for_each_simplex(fun); + + // We could avoid computing maxf, but when !exact rounding errors may cause + // the filtration values to be non-monotonous, so we would need to call + // if (!exact) complex.make_filtration_non_decreasing(); + // which is way more costly than computing maxf. The exact case is already so + // costly that it isn't worth maintaining code without maxf just for it. + // Cech_complex has "free" access to the max of the faces, because + // expansion_with_blockers computes it before the callback. + + // TODO: use a map if complex does not provide key? +} +} // namespace Gudhi::cech_complex + +#endif // MEB_FILTRATION_H_ diff --git a/src/Cech_complex/test/test_cech_complex.cpp b/src/Cech_complex/test/test_cech_complex.cpp index f5980e6d1b..0c36448260 100644 --- a/src/Cech_complex/test/test_cech_complex.cpp +++ b/src/Cech_complex/test/test_cech_complex.cpp @@ -19,6 +19,7 @@ #include // std::max #include +#include // to construct Cech_complex from a OFF file of points #include #include @@ -167,6 +168,13 @@ BOOST_AUTO_TEST_CASE(Cech_complex_for_documentation) { BOOST_CHECK((st2.find({6, 7, 8}) == st2.null_simplex())); BOOST_CHECK((st2.find({3, 5, 7}) == st2.null_simplex())); + + auto st2_save = st2; + st2.reset_filtration(-1); // unnecessary, but ensures we don't cheat + Gudhi::cech_complex::assign_MEB_filtration(Kernel(), st2, points); + for (auto sh : st2.complex_simplex_range()) + st2.assign_filtration(sh, std::sqrt(st2.filtration(sh))); + BOOST_CHECK(st2 == st2_save); // Should only be an approximate test } BOOST_AUTO_TEST_CASE(Cech_complex_from_points) { diff --git a/src/Collapse/include/gudhi/Flag_complex_edge_collapser.h b/src/Collapse/include/gudhi/Flag_complex_edge_collapser.h index d0b3fe4aa7..71679dcf04 100644 --- a/src/Collapse/include/gudhi/Flag_complex_edge_collapser.h +++ b/src/Collapse/include/gudhi/Flag_complex_edge_collapser.h @@ -283,7 +283,7 @@ struct Flag_complex_edge_collapser { }; template R to_range(R&& r) { return std::move(r); } -template R to_range(T&& t) { R r; r.insert(r.end(), t.begin(), t.end()); return r; } +template R to_range(T const& t) { R r; r.insert(r.end(), t.begin(), t.end()); return r; } template auto flag_complex_collapse_edges(FilteredEdgeRange&& edges, Delay&&delay) { diff --git a/src/Subsampling/benchmark/choose_n_farthest_points.cpp b/src/Subsampling/benchmark/choose_n_farthest_points.cpp index c9c453284c..262f8d6f8b 100644 --- a/src/Subsampling/benchmark/choose_n_farthest_points.cpp +++ b/src/Subsampling/benchmark/choose_n_farthest_points.cpp @@ -76,7 +76,16 @@ int main(int argc, char**argv) { << " vs metric " << std::chrono::duration_cast((time_stop2 - time_start2)).count() << " (Boost version " << BOOST_VERSION << ")\n"; if (dists != dists2 || results != results2) { + // Note that with many points, it often happens that 2 points with the same distance are swapped in the output. std::cerr << "Results differ\n"; +#ifdef LOG_DIFF + std::ofstream log_gen("log_gen"); + std::ofstream log_met("log_met"); + for(std::size_t i = 0; i < results.size(); ++i){ + log_gen << dists2[i] << '\t' << results2[i] << '\n'; + log_met << dists [i] << '\t' << results [i] << '\n'; + } +#endif return -1; } #endif diff --git a/src/Subsampling/include/gudhi/choose_n_farthest_points.h b/src/Subsampling/include/gudhi/choose_n_farthest_points.h index 8b57a23c70..eaf91e676d 100644 --- a/src/Subsampling/include/gudhi/choose_n_farthest_points.h +++ b/src/Subsampling/include/gudhi/choose_n_farthest_points.h @@ -320,7 +320,7 @@ void choose_n_farthest_points_metric(Distance dist_, auto handle_neighbor_neighbors = [&](std::size_t ngb) { auto& ngb_info = landmarks[ngb]; - std::remove_if(ngb_info.neighbors.begin(), ngb_info.neighbors.end(), [&](auto near_){ + auto it = std::remove_if(ngb_info.neighbors.begin(), ngb_info.neighbors.end(), [&](auto near_){ std::size_t near = near_.first; FT d = near_.second; // Conservative 3 * radius: we could use the old radii of ngb and near, but not the new ones. @@ -328,6 +328,7 @@ void choose_n_farthest_points_metric(Distance dist_, // Here it is safe to use the new radii. return d >= max_dist(ngb_info.radius, landmarks[near].radius); }); + ngb_info.neighbors.erase(it, ngb_info.neighbors.end()); }; // First update the Voronoi diagram, so we can compute all the updated // radii before pruning neighbor lists. The main drawback is that we have diff --git a/src/cmake/modules/GUDHI_submodules.cmake b/src/cmake/modules/GUDHI_submodules.cmake index 9ede852dc7..9e8124bfc9 100644 --- a/src/cmake/modules/GUDHI_submodules.cmake +++ b/src/cmake/modules/GUDHI_submodules.cmake @@ -3,3 +3,9 @@ set(HERA_INTERNAL_INCLUDE_DIR ${CMAKE_SOURCE_DIR}/ext/hera/include) set(HERA_INCLUDE_DIR ${HERA_INTERNAL_INCLUDE_DIR} CACHE PATH "Directory where one can find hera/{wasserstein.h,bottleneck.h}") # since everything is cleanly under include/hera/, there is no harm always including it include_directories(${HERA_INCLUDE_DIR}) + +if (NOT EXISTS ${HERA_INCLUDE_DIR}/hera/wasserstein.h OR NOT EXISTS ${HERA_INCLUDE_DIR}/hera/bottleneck.h) + message(WARNING "${HERA_INCLUDE_DIR}/hera/{wasserstein.h,bottleneck.h} are not found.\n\ + GUDHI requires this submodules, please consider to launch `git submodule update --init`.\n\ + If hera was installed in a specific directory, you can also consider to specify it to the cmake command with `cmake -DHERA_INCLUDE_DIR=... ...`") +endif() \ No newline at end of file diff --git a/src/cmake/modules/GUDHI_third_party_libraries.cmake b/src/cmake/modules/GUDHI_third_party_libraries.cmake index 4c50fbf727..be89eee80c 100644 --- a/src/cmake/modules/GUDHI_third_party_libraries.cmake +++ b/src/cmake/modules/GUDHI_third_party_libraries.cmake @@ -191,6 +191,7 @@ if (WITH_GUDHI_PYTHON) find_python_module("tensorflow") find_python_module("sphinx_paramlinks") find_python_module("pydata_sphinx_theme") + find_python_module_no_version("sphinxcontrib.bibtex") find_python_module("networkx") endif() diff --git a/src/python/CMakeLists.txt b/src/python/CMakeLists.txt index a7eed460ad..ac48f7eb39 100644 --- a/src/python/CMakeLists.txt +++ b/src/python/CMakeLists.txt @@ -56,6 +56,7 @@ if(PYTHONINTERP_FOUND) set(GUDHI_PYTHON_MODULES "${GUDHI_PYTHON_MODULES}'off_utils', ") set(GUDHI_PYTHON_MODULES "${GUDHI_PYTHON_MODULES}'simplex_tree', ") set(GUDHI_PYTHON_MODULES "${GUDHI_PYTHON_MODULES}'simplex_tree_multi', ") + set(GUDHI_PYTHON_MODULES "${GUDHI_PYTHON_MODULES}'edge_collapse', ") set(GUDHI_PYTHON_MODULES "${GUDHI_PYTHON_MODULES}'rips_complex', ") set(GUDHI_PYTHON_MODULES "${GUDHI_PYTHON_MODULES}'cubical_complex', ") set(GUDHI_PYTHON_MODULES "${GUDHI_PYTHON_MODULES}'periodic_cubical_complex', ") @@ -122,9 +123,12 @@ if(PYTHONINTERP_FOUND) add_gudhi_debug_info("Sphinx-paramlinks version ${SPHINX_PARAMLINKS_VERSION}") endif() if(PYDATA_SPHINX_THEME_FOUND) - # Does not have a version number... add_gudhi_debug_info("pydata_sphinx_theme version ${PYDATA_SPHINX_THEME_VERSION}") endif() + if(SPHINXCONTRIB.BIBTEX_FOUND) + # Does not have a version number... + add_gudhi_debug_info("sphinxcontrib-bibtex found") + endif() if(NETWORKX_FOUND) add_gudhi_debug_info("NetworkX version ${NETWORKX_VERSION}") endif() @@ -187,6 +191,7 @@ if(PYTHONINTERP_FOUND) set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'alpha_complex', ") endif () set(GUDHI_PYBIND11_MODULES "${GUDHI_PYBIND11_MODULES}'_pers_cub_low_dim', ") + set(GUDHI_PYBIND11_MODULES "${GUDHI_PYBIND11_MODULES}'_edge_collapse', ") # from windows vcpkg eigen 3.4.0#2 : build fails with # error C2440: '': cannot convert from 'Eigen::EigenBase::Index' to '__gmp_expr' @@ -325,6 +330,7 @@ if(PYTHONINTERP_FOUND) file(COPY "gudhi/hera/__init__.py" DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/gudhi/hera") file(COPY "gudhi/datasets" DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/gudhi" FILES_MATCHING PATTERN "*.py") file(COPY "gudhi/sklearn" DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/gudhi/") + file(COPY "gudhi/edge_collapse.py" DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/gudhi") # Some files for pip package @@ -348,66 +354,71 @@ if(PYTHONINTERP_FOUND) # Make it first as sphinx test is by far the longest test which is nice when testing in parallel if(SPHINX_PATH) if(SPHINX_PARAMLINKS_FOUND) - if(PYDATA_SPHINX_THEME_FOUND) - if(MATPLOTLIB_FOUND) - if(NUMPY_FOUND) - if(SCIPY_FOUND) - if(SKLEARN_FOUND) - if(OT_FOUND) - if(NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 5.1.0) - set (GUDHI_SPHINX_MESSAGE "Generating API documentation with Sphinx in ${CMAKE_CURRENT_BINARY_DIR}/sphinx/") - # User warning - Sphinx is a static pages generator, and configured to work fine with user_version - # Images and biblio warnings because not found on developer version - if (GUDHI_PYTHON_PATH STREQUAL "src/python") - set (GUDHI_SPHINX_MESSAGE "${GUDHI_SPHINX_MESSAGE} \n WARNING : Sphinx is configured for user version, you run it on developer version. Images and biblio will miss") - endif() - # sphinx target requires gudhi.so, because conf.py reads gudhi version from it - add_custom_target(sphinx - WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}/doc - COMMAND ${CMAKE_COMMAND} -E env "${GUDHI_PYTHON_PATH_ENV}" - ${SPHINX_PATH} -b html ${CMAKE_CURRENT_SOURCE_DIR}/doc ${CMAKE_CURRENT_BINARY_DIR}/sphinx - DEPENDS "${CMAKE_CURRENT_BINARY_DIR}/gudhi.so" - COMMENT "${GUDHI_SPHINX_MESSAGE}" VERBATIM) - add_test(NAME sphinx_py_test - WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} - COMMAND ${CMAKE_COMMAND} -E env "${GUDHI_PYTHON_PATH_ENV}" - ${SPHINX_PATH} -b doctest ${CMAKE_CURRENT_SOURCE_DIR}/doc ${CMAKE_CURRENT_BINARY_DIR}/doctest) - # Set missing or not modules - set(GUDHI_MODULES ${GUDHI_MODULES} "python-documentation" CACHE INTERNAL "GUDHI_MODULES") - else(NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 5.1.0) - message("++ Python documentation module will not be compiled because it requires a Eigen3 and CGAL version >= 5.1.0") + if(SPHINXCONTRIB.BIBTEX_FOUND) + if(PYDATA_SPHINX_THEME_FOUND) + if(MATPLOTLIB_FOUND) + if(NUMPY_FOUND) + if(SCIPY_FOUND) + if(SKLEARN_FOUND) + if(OT_FOUND) + if(NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 5.1.0) + set (GUDHI_SPHINX_MESSAGE "Generating API documentation with Sphinx in ${CMAKE_CURRENT_BINARY_DIR}/sphinx/") + # User warning - Sphinx is a static pages generator, and configured to work fine with user_version + # Images and biblio warnings because not found on developer version + if (GUDHI_PYTHON_PATH STREQUAL "src/python") + set (GUDHI_SPHINX_MESSAGE "${GUDHI_SPHINX_MESSAGE} \n WARNING : Sphinx is configured for user version, you run it on developer version. Images and biblio will miss") + endif() + # sphinx target requires gudhi.so, because conf.py reads gudhi version from it + add_custom_target(sphinx + WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}/doc + COMMAND ${CMAKE_COMMAND} -E env "${GUDHI_PYTHON_PATH_ENV}" + ${SPHINX_PATH} -b html ${CMAKE_CURRENT_SOURCE_DIR}/doc ${CMAKE_CURRENT_BINARY_DIR}/sphinx + DEPENDS "${CMAKE_CURRENT_BINARY_DIR}/gudhi.so" + COMMENT "${GUDHI_SPHINX_MESSAGE}" VERBATIM) + add_test(NAME sphinx_py_test + WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} + COMMAND ${CMAKE_COMMAND} -E env "${GUDHI_PYTHON_PATH_ENV}" + ${SPHINX_PATH} -b doctest ${CMAKE_CURRENT_SOURCE_DIR}/doc ${CMAKE_CURRENT_BINARY_DIR}/doctest) + # Set missing or not modules + set(GUDHI_MODULES ${GUDHI_MODULES} "python-documentation" CACHE INTERNAL "GUDHI_MODULES") + else(NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 5.1.0) + message("++ Python documentation module will not be compiled because it requires a Eigen3 and CGAL version >= 5.1.0") + set(GUDHI_MISSING_MODULES ${GUDHI_MISSING_MODULES} "python-documentation" CACHE INTERNAL "GUDHI_MISSING_MODULES") + endif(NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 5.1.0) + else(OT_FOUND) + message("++ Python documentation module will not be compiled because POT was not found") set(GUDHI_MISSING_MODULES ${GUDHI_MISSING_MODULES} "python-documentation" CACHE INTERNAL "GUDHI_MISSING_MODULES") - endif(NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 5.1.0) - else(OT_FOUND) - message("++ Python documentation module will not be compiled because POT was not found") + endif(OT_FOUND) + else(SKLEARN_FOUND) + message("++ Python documentation module will not be compiled because scikit-learn was not found") set(GUDHI_MISSING_MODULES ${GUDHI_MISSING_MODULES} "python-documentation" CACHE INTERNAL "GUDHI_MISSING_MODULES") - endif(OT_FOUND) - else(SKLEARN_FOUND) - message("++ Python documentation module will not be compiled because scikit-learn was not found") + endif(SKLEARN_FOUND) + else(SCIPY_FOUND) + message("++ Python documentation module will not be compiled because scipy was not found") set(GUDHI_MISSING_MODULES ${GUDHI_MISSING_MODULES} "python-documentation" CACHE INTERNAL "GUDHI_MISSING_MODULES") - endif(SKLEARN_FOUND) - else(SCIPY_FOUND) - message("++ Python documentation module will not be compiled because scipy was not found") + endif(SCIPY_FOUND) + else(NUMPY_FOUND) + message("++ Python documentation module will not be compiled because numpy was not found") set(GUDHI_MISSING_MODULES ${GUDHI_MISSING_MODULES} "python-documentation" CACHE INTERNAL "GUDHI_MISSING_MODULES") - endif(SCIPY_FOUND) - else(NUMPY_FOUND) - message("++ Python documentation module will not be compiled because numpy was not found") + endif(NUMPY_FOUND) + else(MATPLOTLIB_FOUND) + message("++ Python documentation module will not be compiled because matplotlib was not found") set(GUDHI_MISSING_MODULES ${GUDHI_MISSING_MODULES} "python-documentation" CACHE INTERNAL "GUDHI_MISSING_MODULES") - endif(NUMPY_FOUND) - else(MATPLOTLIB_FOUND) - message("++ Python documentation module will not be compiled because matplotlib was not found") + endif(MATPLOTLIB_FOUND) + else(PYDATA_SPHINX_THEME_FOUND) + message("++ Python documentation module will not be compiled because pydata_sphinx_theme was not found") set(GUDHI_MISSING_MODULES ${GUDHI_MISSING_MODULES} "python-documentation" CACHE INTERNAL "GUDHI_MISSING_MODULES") - endif(MATPLOTLIB_FOUND) - else(PYDATA_SPHINX_THEME_FOUND) - message("++ Python documentation module will not be compiled because pydata_sphinx_theme was not found") + endif(PYDATA_SPHINX_THEME_FOUND) + else(SPHINXCONTRIB.BIBTEX_FOUND) + message("++ Python documentation module will not be compiled because sphinxcontrib-bibtex was not found") set(GUDHI_MISSING_MODULES ${GUDHI_MISSING_MODULES} "python-documentation" CACHE INTERNAL "GUDHI_MISSING_MODULES") - endif(PYDATA_SPHINX_THEME_FOUND) + endif(SPHINXCONTRIB.BIBTEX_FOUND) else(SPHINX_PARAMLINKS_FOUND) message("++ Python documentation module will not be compiled because sphinxcontrib-paramlinks was not found") set(GUDHI_MISSING_MODULES ${GUDHI_MISSING_MODULES} "python-documentation" CACHE INTERNAL "GUDHI_MISSING_MODULES") endif(SPHINX_PARAMLINKS_FOUND) else(SPHINX_PATH) - message("++ Python documentation module will not be compiled because sphinx and sphinxcontrib-bibtex were not found") + message("++ Python documentation module will not be compiled because sphinx was not found") set(GUDHI_MISSING_MODULES ${GUDHI_MISSING_MODULES} "python-documentation" CACHE INTERNAL "GUDHI_MISSING_MODULES") endif(SPHINX_PATH) @@ -565,6 +576,11 @@ if(PYTHONINTERP_FOUND) add_gudhi_py_test(test_simplex_generators) add_gudhi_py_test(test_simplex_tree_serialization) + # Edge collapse + if(SCIPY_FOUND) + add_gudhi_py_test(test_collapse_edges) + endif() + # Subsampling add_gudhi_py_test(test_subsampling) diff --git a/src/python/doc/edge_collapse.rst b/src/python/doc/edge_collapse.rst new file mode 100644 index 0000000000..9d727b7aa6 --- /dev/null +++ b/src/python/doc/edge_collapse.rst @@ -0,0 +1,8 @@ +:orphan: + +.. To get rid of WARNING: document isn't included in any toctree + +Edge collapse +============= + +.. autofunction:: gudhi.collapse_edges diff --git a/src/python/doc/edge_collapse_sum.inc b/src/python/doc/edge_collapse_sum.inc new file mode 100644 index 0000000000..bc6c6d3acd --- /dev/null +++ b/src/python/doc/edge_collapse_sum.inc @@ -0,0 +1,18 @@ +.. table:: + :widths: 30 40 30 + + +-----------------------------------------------------------------+---------------------------------------------------------------------+-----------------------------------------------+ + | .. figure:: | Edge collapse is able to reduce any flag filtration to a smaller | :Author: Siddharth Pritam, Marc Glisse | + | ../../doc/Collapse/dominated_edge.png | flag filtration with the same persistence, using only the | | + | :figclass: align-center | 1-skeleton of a simplicial complex. | :Since: GUDHI 3.9.0 | + | | The reduction is exact and the persistent homology of the reduced | | + | | sequence is identical to the persistent homology of the input | | + | | sequence. The resulting method is simple and extremely efficient. | | + | | | | + | | Computation of edge collapse and persistent homology of a filtered | | + | | flag complex via edge collapse as described in | | + | | :cite:`edgecollapsearxiv`. | | + | | | | + +-----------------------------------------------------------------+---------------------------------------------------------------------+-----------------------------------------------+ + | | * :doc:`edge_collapse` | + +-----------------------------------------------------------------+---------------------------------------------------------------------------------------------------------------------+ diff --git a/src/python/doc/index.rst b/src/python/doc/index.rst index 35f4ba466a..7aa7b52527 100644 --- a/src/python/doc/index.rst +++ b/src/python/doc/index.rst @@ -35,6 +35,11 @@ Rips complex .. include:: rips_complex_sum.inc +Edge collapse +============= + +.. include:: edge_collapse_sum.inc + Witness complex =============== diff --git a/src/python/gudhi/_edge_collapse.cc b/src/python/gudhi/_edge_collapse.cc new file mode 100644 index 0000000000..40e176f76a --- /dev/null +++ b/src/python/gudhi/_edge_collapse.cc @@ -0,0 +1,67 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): Marc Glisse + * + * Copyright (C) 2023 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#include +#include +#include + +#include +#include +#include + +#include + +namespace py = pybind11; + +template +py::object collapse(py::array_t is, py::array_t js, py::array_t fs, int nb_iterations) { + typedef std::tuple Filtered_edge; + typedef std::vector Edges; + py::buffer_info bufi = is.request(); + py::buffer_info bufj = js.request(); + py::buffer_info buff = fs.request(); + if (bufi.ndim != 1 || bufj.ndim != 1 || buff.ndim != 1) + throw std::runtime_error("Input must be 1-dimensional arrays"); + if (bufi.shape[0] != bufj.shape[0] || bufi.shape[0] != buff.shape[0]) + throw std::runtime_error("Input arrays must have the same size"); + if (buff.shape[0] == 0) { + py::array_t indices({{ 2, 0 }}, {{ 0, 0 }}); + py::array_t filtrs; + return py::make_tuple(std::move(indices), std::move(filtrs)); + } + auto& edges = *new Edges(); + { + py::gil_scoped_release release; + Index n_edges = static_cast(bufi.shape[0]); + edges.reserve(n_edges); + auto strides_i = bufi.strides[0]; + auto strides_j = bufj.strides[0]; + auto strides_f = buff.strides[0]; + for (Index k = 0; k < n_edges; ++k) { + Index i = *reinterpret_cast(static_cast(bufi.ptr) + k * strides_i); + Index j = *reinterpret_cast(static_cast(bufj.ptr) + k * strides_j); + Filtr f = *reinterpret_cast(static_cast(buff.ptr) + k * strides_f); + edges.emplace_back(i, j, f); + } + for (int k = 0; k < nb_iterations; ++k) { + edges = Gudhi::collapse::flag_complex_collapse_edges(std::move(edges), [](auto const&d){return d;}); + } + } + py::capsule owner(&edges, [](void*p){ delete reinterpret_cast(p); }); + const auto offset = reinterpret_cast(&std::get<1>(edges[0])) - reinterpret_cast(&std::get<0>(edges[0])); + py::array_t indices({{ 2, static_cast(edges.size()) }}, {{ offset, sizeof(Filtered_edge) }}, &std::get<0>(edges[0]), owner); + py::array_t filtrs ({{ static_cast(edges.size()) }}, {{ sizeof(Filtered_edge) }}, &std::get<2>(edges[0]), owner); + return py::make_tuple(std::move(indices), std::move(filtrs)); +} + +PYBIND11_MODULE(_edge_collapse, m) { + m.def("_collapse_edges", collapse, py::arg("i").noconvert(), py::arg("j").noconvert(), py::arg("f").noconvert(), py::arg("nb_iterations")=1); + m.def("_collapse_edges", collapse, py::arg("i"), py::arg("j"), py::arg("f"), py::arg("nb_iterations")=1); +} diff --git a/src/python/gudhi/edge_collapse.py b/src/python/gudhi/edge_collapse.py new file mode 100644 index 0000000000..e1ac90c1c6 --- /dev/null +++ b/src/python/gudhi/edge_collapse.py @@ -0,0 +1,35 @@ +# This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. +# See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. +# Author(s): Marc Glisse +# +# Copyright (C) 2023 Inria +# +# Modification(s): +# - YYYY/MM Author: Description of the modification + +from ._edge_collapse import _collapse_edges + + +def collapse_edges(input_edges, nb_iterations=1): + """Simplify a clique filtration, preserving its persistent homology. + The clique filtration is represented as a sparse weighted graph, giving for each edge its 2 vertices + and its filtration value, and the output is in the same format. An edge may be represented arbitrarily as (i, j) + or (j, i). Listing the same edge twice, or a self-loop, is undefined. + The cliques of the graph composed of the edges with filtration value less than some arbitrary value implicitly + define a flag complex, so the weighted graph describes a flag filtration. + This function outputs another flag filtration which is guaranteed to have the same persistent homology as the input. + The algorithm works by collapsing edges, as described in :cite:`edgecollapsearxiv`. + Note that this simplification is independent of the filtration values of the vertices. + The output has the same shape as the input, which is presumed to be (N, N) where all vertices have index + less than N, since the simplification does not affect vertices. + + :param input_edges: Input weighted graph. + :type input_edges: scipy.sparse.coo_matrix + :param nb_iterations: The number of times we apply the algorithm. Default is 1. + :type nb_iterations: int + :rtype: scipy.sparse.coo_matrix + """ + from scipy.sparse import coo_matrix + + r = _collapse_edges(input_edges.row, input_edges.col, input_edges.data, nb_iterations) + return coo_matrix((r[1], (r[0][0], r[0][1])), input_edges.shape) diff --git a/src/python/test/test_collapse_edges.py b/src/python/test/test_collapse_edges.py new file mode 100644 index 0000000000..67e31acbf0 --- /dev/null +++ b/src/python/test/test_collapse_edges.py @@ -0,0 +1,38 @@ +""" This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + Author(s): Marc Glisse + + Copyright (C) 2023 Inria + + Modification(s): + - YYYY/MM Author: Description of the modification +""" + +from gudhi import collapse_edges +import numpy as np +from scipy.sparse import coo_matrix +import pytest + + +def test_collapse_empty(): + x = coo_matrix(([], ([], [])), (3, 3)) + xo = collapse_edges(x) + assert xo.shape == (3, 3) + assert len(xo.data) == 0 + + +def test_collapse(): + x = coo_matrix(([0.1, 0.2, 0.3, 0.4, 0.6], ([0, 1, 2, 3, 1], [1, 2, 3, 0, 3])), (6, 6)) + xo = collapse_edges(x) + assert xo.shape == (6, 6) + assert len(xo.data) == 5 + + x = coo_matrix(([0.1, 0.2, 0.3, 0.4, -0.6], ([0, 1, 2, 3, 1], [1, 2, 3, 0, 3])), (4, 4)) + xo = collapse_edges(x) + assert xo.shape == (4, 4) + assert len(xo.data) == 3 + assert xo.dtype == np.dtype("float64") + + x = coo_matrix((np.array([0.1], dtype="float32"), ([0], [1])), (2, 2)) + xo = collapse_edges(x) + assert xo.dtype == np.dtype("float32")