From 4319c65cc621cab612f7e3732c5266c2cb552846 Mon Sep 17 00:00:00 2001 From: Marc Glisse Date: Sun, 2 Jul 2023 20:59:36 +0200 Subject: [PATCH 01/17] python wrapper for edge collapse --- .../gudhi/Flag_complex_edge_collapser.h | 2 +- src/python/CMakeLists.txt | 1 + src/python/gudhi/_edge_collapse.cc | 59 +++++++++++++++++++ 3 files changed, 61 insertions(+), 1 deletion(-) create mode 100644 src/python/gudhi/_edge_collapse.cc 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/python/CMakeLists.txt b/src/python/CMakeLists.txt index 86a409b6b8..e79c40be7f 100644 --- a/src/python/CMakeLists.txt +++ b/src/python/CMakeLists.txt @@ -185,6 +185,7 @@ if(PYTHONINTERP_FOUND) set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'alpha_complex', ") endif () set(GUDHI_PYBIND11_MODULES "${GUDHI_PYBIND11_MODULES}'_persline', ") + 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' diff --git a/src/python/gudhi/_edge_collapse.cc b/src/python/gudhi/_edge_collapse.cc new file mode 100644 index 0000000000..7e0ccd3a32 --- /dev/null +++ b/src/python/gudhi/_edge_collapse.cc @@ -0,0 +1,59 @@ +/* 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 + +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"); + auto& edges = *new Edges();; + edges.reserve(bufi.shape[0]); + Index n_edges = static_cast(bufi.shape[0]); + 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;}); + } + auto res = py::array(py::cast(std::move(edges))); + 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(indices, 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); +} From 2456ab4e8273a663f63b8771b188de1219995d51 Mon Sep 17 00:00:00 2001 From: Marc Glisse Date: Sun, 2 Jul 2023 23:15:18 +0200 Subject: [PATCH 02/17] collapse_edges python API --- src/python/CMakeLists.txt | 5 +++ src/python/doc/edge_collapse.rst | 8 +++++ src/python/doc/persistent_cohomology_sum.inc | 1 + src/python/gudhi/_edge_collapse.cc | 8 +++-- src/python/gudhi/edge_collapse.py | 35 ++++++++++++++++++ src/python/test/test_collapse_edges.py | 38 ++++++++++++++++++++ 6 files changed, 93 insertions(+), 2 deletions(-) create mode 100644 src/python/doc/edge_collapse.rst create mode 100644 src/python/gudhi/edge_collapse.py create mode 100644 src/python/test/test_collapse_edges.py diff --git a/src/python/CMakeLists.txt b/src/python/CMakeLists.txt index e79c40be7f..75e0cf121c 100644 --- a/src/python/CMakeLists.txt +++ b/src/python/CMakeLists.txt @@ -55,6 +55,7 @@ if(PYTHONINTERP_FOUND) # Cython modules 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}'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', ") @@ -320,6 +321,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 @@ -560,6 +562,9 @@ if(PYTHONINTERP_FOUND) add_gudhi_py_test(test_simplex_generators) add_gudhi_py_test(test_simplex_tree_serialization) + # Edge collapse + add_gudhi_py_test(test_collapse_edges) + # 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/persistent_cohomology_sum.inc b/src/python/doc/persistent_cohomology_sum.inc index 58e44b8a96..fb2d99ecf4 100644 --- a/src/python/doc/persistent_cohomology_sum.inc +++ b/src/python/doc/persistent_cohomology_sum.inc @@ -21,4 +21,5 @@ | | * :doc:`simplex_tree_ref` | | | * :doc:`cubical_complex_ref` | | | * :doc:`periodic_cubical_complex_ref` | + | | * :doc:`edge_collapse` | +-----------------------------------------------------------------+-----------------------------------------------------------------------------------------------------------------------+ diff --git a/src/python/gudhi/_edge_collapse.cc b/src/python/gudhi/_edge_collapse.cc index 7e0ccd3a32..5d500e413d 100644 --- a/src/python/gudhi/_edge_collapse.cc +++ b/src/python/gudhi/_edge_collapse.cc @@ -30,7 +30,12 @@ py::object collapse(py::array_t is, py::array_t js, 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(); edges.reserve(bufi.shape[0]); Index n_edges = static_cast(bufi.shape[0]); auto strides_i = bufi.strides[0]; @@ -45,7 +50,6 @@ py::object collapse(py::array_t is, py::array_t js, py::array_t(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); diff --git a/src/python/gudhi/edge_collapse.py b/src/python/gudhi/edge_collapse.py new file mode 100644 index 0000000000..c0a0b8271d --- /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): Vincent Rouvreau +# +# 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, perserving 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) + 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") From ce5946dfacda8e91902aec57b3a151846e6e6a43 Mon Sep 17 00:00:00 2001 From: Marc Glisse Date: Mon, 3 Jul 2023 10:09:45 +0200 Subject: [PATCH 03/17] GIL --- src/python/gudhi/_edge_collapse.cc | 31 ++++++++++++++++-------------- 1 file changed, 17 insertions(+), 14 deletions(-) diff --git a/src/python/gudhi/_edge_collapse.cc b/src/python/gudhi/_edge_collapse.cc index 5d500e413d..3a9edc6673 100644 --- a/src/python/gudhi/_edge_collapse.cc +++ b/src/python/gudhi/_edge_collapse.cc @@ -36,25 +36,28 @@ py::object collapse(py::array_t is, py::array_t js, py::array_t(bufi.shape[0]); - 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::gil_scoped_release release; + edges.reserve(bufi.shape[0]); + Index n_edges = static_cast(bufi.shape[0]); + 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(indices, filtrs); + return py::make_tuple(std::move(indices), std::move(filtrs)); } PYBIND11_MODULE(_edge_collapse, m) { From a99da7cb89a947e5da6e1af9a1fd22880422a7e7 Mon Sep 17 00:00:00 2001 From: Marc Glisse Date: Mon, 28 Aug 2023 11:26:12 +0200 Subject: [PATCH 04/17] Review comments --- src/python/CMakeLists.txt | 4 +++- src/python/gudhi/_edge_collapse.cc | 3 ++- src/python/gudhi/edge_collapse.py | 4 ++-- 3 files changed, 7 insertions(+), 4 deletions(-) diff --git a/src/python/CMakeLists.txt b/src/python/CMakeLists.txt index 75e0cf121c..7f7c63e04d 100644 --- a/src/python/CMakeLists.txt +++ b/src/python/CMakeLists.txt @@ -563,7 +563,9 @@ if(PYTHONINTERP_FOUND) add_gudhi_py_test(test_simplex_tree_serialization) # Edge collapse - add_gudhi_py_test(test_collapse_edges) + if(SCIPY_FOUND) + add_gudhi_py_test(test_collapse_edges) + endif() # Subsampling add_gudhi_py_test(test_subsampling) diff --git a/src/python/gudhi/_edge_collapse.cc b/src/python/gudhi/_edge_collapse.cc index 3a9edc6673..40e176f76a 100644 --- a/src/python/gudhi/_edge_collapse.cc +++ b/src/python/gudhi/_edge_collapse.cc @@ -10,6 +10,7 @@ #include #include +#include #include #include @@ -38,8 +39,8 @@ py::object collapse(py::array_t is, py::array_t js, py::array_t(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]; diff --git a/src/python/gudhi/edge_collapse.py b/src/python/gudhi/edge_collapse.py index c0a0b8271d..ba9edad6ea 100644 --- a/src/python/gudhi/edge_collapse.py +++ b/src/python/gudhi/edge_collapse.py @@ -1,6 +1,6 @@ # This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. # See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. -# Author(s): Vincent Rouvreau +# Author(s): Marc Glisse # # Copyright (C) 2023 Inria # @@ -11,7 +11,7 @@ def collapse_edges(input_edges, nb_iterations=1): - """Simplify a clique filtration, perserving its persistent homology. + """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. From a04a8fd0076eb5706f87dab4eb8622295f399058 Mon Sep 17 00:00:00 2001 From: Marc Glisse Date: Tue, 29 Aug 2023 18:18:14 +0200 Subject: [PATCH 05/17] missing argument --- src/python/gudhi/edge_collapse.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/python/gudhi/edge_collapse.py b/src/python/gudhi/edge_collapse.py index ba9edad6ea..e1ac90c1c6 100644 --- a/src/python/gudhi/edge_collapse.py +++ b/src/python/gudhi/edge_collapse.py @@ -31,5 +31,5 @@ def collapse_edges(input_edges, nb_iterations=1): """ from scipy.sparse import coo_matrix - r = _collapse_edges(input_edges.row, input_edges.col, input_edges.data) + 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) From da633e6ac06868e9031954e9c5bc52be12166da7 Mon Sep 17 00:00:00 2001 From: Marc Glisse Date: Sun, 23 Jul 2023 13:27:47 +0200 Subject: [PATCH 06/17] New assign_MEB_filtration --- .../include/gudhi/Cech_complex_blocker.h | 2 +- .../include/gudhi/MEB_filtration.h | 95 +++++++++++++++++++ src/Cech_complex/test/test_cech_complex.cpp | 8 ++ 3 files changed, 104 insertions(+), 1 deletion(-) create mode 100644 src/Cech_complex/include/gudhi/MEB_filtration.h 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..ba8e8f2e87 --- /dev/null +++ b/src/Cech_complex/include/gudhi/MEB_filtration.h @@ -0,0 +1,95 @@ +/* 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 + * + * 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. + * + * Applied on a Cech complex, it recomputes the same values (squared). Applied on a Delaunay triangulation, it computes the Delaunay-Cech filtration. + * + * \tparam Kernel CGAL kernel: either Epick_d or Epeck_d. + * \tparam SimplicialComplexForMEB Simplex_tree with Simplex_key an integer type large enough to store the number of simplices. + * \tparam PointRange Random access range of `Kernel::Point_d`. + * + * @param[in] points Embedding of the vertices of the complex. + * @param[in] complex The simplicial 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) { + // 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()); // TODO: use a map if complex does not provide key? + complex.assign_filtration(sh, std::max(cvt(r), Filtration_value(0))); + cache_.emplace_back(std::move(m), std::move(r)); + } else { + for (auto face_opposite_vertex : complex.boundary_opposite_vertex_simplex_range(sh)) { + 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; + complex.assign_key(sh, key); + complex.assign_filtration(sh, complex.filtration(face_opposite_vertex.first)); + return; + } + // 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); + complex.assign_key(sh, cache_.size()); // TODO: use a map if complex does not provide key? + complex.assign_filtration(sh, cvt(r)); + cache_.emplace_back(std::move(c), std::move(r)); + } + }; + complex.for_each_simplex_with_dim(fun); + // TODO: when !exact, how do we ensure that the value is indeed larger + // than for the faces? Cech_complex also looks at the max of the faces, + // given by expansion. We would need to do compute this max as well, and + // in particular not short-circuit when we find a good face. For now, the + // simplest: + if (!exact) complex.make_filtration_non_decreasing(); +} +} // 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) { From aea8f4a68780ade1fc504614bd118ebd6cd11542 Mon Sep 17 00:00:00 2001 From: Marc Glisse Date: Thu, 14 Sep 2023 17:51:50 +0200 Subject: [PATCH 07/17] Finish renaming for_each_simplex --- src/Cech_complex/include/gudhi/MEB_filtration.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Cech_complex/include/gudhi/MEB_filtration.h b/src/Cech_complex/include/gudhi/MEB_filtration.h index ba8e8f2e87..1f02750a5f 100644 --- a/src/Cech_complex/include/gudhi/MEB_filtration.h +++ b/src/Cech_complex/include/gudhi/MEB_filtration.h @@ -82,7 +82,7 @@ void assign_MEB_filtration(Kernel&&k, SimplicialComplexForMEB& complex, PointRan cache_.emplace_back(std::move(c), std::move(r)); } }; - complex.for_each_simplex_with_dim(fun); + complex.for_each_simplex(fun); // TODO: when !exact, how do we ensure that the value is indeed larger // than for the faces? Cech_complex also looks at the max of the faces, // given by expansion. We would need to do compute this max as well, and From cf61a1f589df75bc534ac3b13b1952d245760249 Mon Sep 17 00:00:00 2001 From: Vincent Rouvreau Date: Fri, 15 Sep 2023 18:04:33 +0200 Subject: [PATCH 08/17] Python documentation requires sphinxcontrib-bibtex to be built --- .../modules/GUDHI_third_party_libraries.cmake | 1 + src/python/CMakeLists.txt | 104 ++++++++++-------- 2 files changed, 57 insertions(+), 48 deletions(-) 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 6540f2a785..7efc9915f5 100644 --- a/src/python/CMakeLists.txt +++ b/src/python/CMakeLists.txt @@ -121,9 +121,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() @@ -346,60 +349,65 @@ 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") From 1488c9bed79a0371cff48f12f51db0ea8a6ac840 Mon Sep 17 00:00:00 2001 From: Marc Glisse Date: Sun, 17 Sep 2023 13:52:13 +0200 Subject: [PATCH 09/17] Avoid make_filtration_non_decreasing --- .../include/gudhi/MEB_filtration.h | 66 ++++++++++++------- 1 file changed, 41 insertions(+), 25 deletions(-) diff --git a/src/Cech_complex/include/gudhi/MEB_filtration.h b/src/Cech_complex/include/gudhi/MEB_filtration.h index 1f02750a5f..d4f9c65c39 100644 --- a/src/Cech_complex/include/gudhi/MEB_filtration.h +++ b/src/Cech_complex/include/gudhi/MEB_filtration.h @@ -48,8 +48,8 @@ void assign_MEB_filtration(Kernel&&k, SimplicialComplexForMEB& complex, PointRan auto fun = [&](Simplex_handle sh, int dim){ if (dim == 0) complex.assign_filtration(sh, 0); else if (dim == 1) { - // Vertex_handle u = sh->first; - // Vertex_handle v = self_siblings(sh)->parent(); + // 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; @@ -58,37 +58,53 @@ void assign_MEB_filtration(Kernel&&k, SimplicialComplexForMEB& complex, PointRan 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()); // TODO: use a map if complex does not provide key? + 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)) { - 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; - complex.assign_key(sh, key); - complex.assign_filtration(sh, complex.filtration(face_opposite_vertex.first)); - return; + 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 + } } - // 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); - complex.assign_key(sh, cache_.size()); // TODO: use a map if complex does not provide key? - complex.assign_filtration(sh, cvt(r)); - cache_.emplace_back(std::move(c), std::move(r)); + 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); - // TODO: when !exact, how do we ensure that the value is indeed larger - // than for the faces? Cech_complex also looks at the max of the faces, - // given by expansion. We would need to do compute this max as well, and - // in particular not short-circuit when we find a good face. For now, the - // simplest: - if (!exact) complex.make_filtration_non_decreasing(); + + // 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 From a12b15e773582c177c96c756ff5d521d48962a77 Mon Sep 17 00:00:00 2001 From: Marc Glisse Date: Sun, 17 Sep 2023 14:35:51 +0200 Subject: [PATCH 10/17] Improve doc --- .../concept/SimplicialComplexForMEB.h | 58 +++++++++++++++++++ src/Cech_complex/doc/Intro_cech_complex.h | 7 ++- .../include/gudhi/MEB_filtration.h | 5 +- 3 files changed, 67 insertions(+), 3 deletions(-) create mode 100644 src/Cech_complex/concept/SimplicialComplexForMEB.h diff --git a/src/Cech_complex/concept/SimplicialComplexForMEB.h b/src/Cech_complex/concept/SimplicialComplexForMEB.h new file mode 100644 index 0000000000..f638867a8c --- /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/MEB_filtration.h b/src/Cech_complex/include/gudhi/MEB_filtration.h index d4f9c65c39..c1eefe8dbb 100644 --- a/src/Cech_complex/include/gudhi/MEB_filtration.h +++ b/src/Cech_complex/include/gudhi/MEB_filtration.h @@ -16,6 +16,7 @@ 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. @@ -23,11 +24,11 @@ namespace Gudhi::cech_complex { * Applied on a Cech complex, it recomputes the same values (squared). Applied on a Delaunay triangulation, it computes the Delaunay-Cech filtration. * * \tparam Kernel CGAL kernel: either Epick_d or Epeck_d. - * \tparam SimplicialComplexForMEB Simplex_tree with Simplex_key an integer type large enough to store the number of simplices. * \tparam PointRange Random access range of `Kernel::Point_d`. * - * @param[in] points Embedding of the vertices of the complex. + * @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. */ From 792db83a04b6948eca06142ce3578802b018aac7 Mon Sep 17 00:00:00 2001 From: Marc Glisse Date: Sun, 17 Sep 2023 14:49:31 +0200 Subject: [PATCH 11/17] Doc tweak --- src/Cech_complex/concept/SimplicialComplexForMEB.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Cech_complex/concept/SimplicialComplexForMEB.h b/src/Cech_complex/concept/SimplicialComplexForMEB.h index f638867a8c..19a0a90818 100644 --- a/src/Cech_complex/concept/SimplicialComplexForMEB.h +++ b/src/Cech_complex/concept/SimplicialComplexForMEB.h @@ -35,7 +35,7 @@ struct SimplicialComplexForMEB { /** \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`. */ + /** \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); From 731d3876b282140a50963d768b2d3c7ee614d9c1 Mon Sep 17 00:00:00 2001 From: Marc Glisse Date: Sun, 17 Sep 2023 15:32:48 +0200 Subject: [PATCH 12/17] Move doc to its own section in index --- src/python/doc/edge_collapse_sum.inc | 18 ++++++++++++++++++ src/python/doc/index.rst | 5 +++++ src/python/doc/persistent_cohomology_sum.inc | 1 - 3 files changed, 23 insertions(+), 1 deletion(-) create mode 100644 src/python/doc/edge_collapse_sum.inc 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/doc/persistent_cohomology_sum.inc b/src/python/doc/persistent_cohomology_sum.inc index fb2d99ecf4..58e44b8a96 100644 --- a/src/python/doc/persistent_cohomology_sum.inc +++ b/src/python/doc/persistent_cohomology_sum.inc @@ -21,5 +21,4 @@ | | * :doc:`simplex_tree_ref` | | | * :doc:`cubical_complex_ref` | | | * :doc:`periodic_cubical_complex_ref` | - | | * :doc:`edge_collapse` | +-----------------------------------------------------------------+-----------------------------------------------------------------------------------------------------------------------+ From 04e5bd3176b09763c4afcb9630c2ab9ca15b16e1 Mon Sep 17 00:00:00 2001 From: Vincent Rouvreau Date: Mon, 18 Sep 2023 14:05:00 +0200 Subject: [PATCH 13/17] code review: old error message was not up to date --- src/python/CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/python/CMakeLists.txt b/src/python/CMakeLists.txt index 7efc9915f5..5bb2f59677 100644 --- a/src/python/CMakeLists.txt +++ b/src/python/CMakeLists.txt @@ -413,7 +413,7 @@ if(PYTHONINTERP_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) From d24894eadc4540b5c6c30261d1c8a48a6802f728 Mon Sep 17 00:00:00 2001 From: Marc Glisse Date: Mon, 18 Sep 2023 20:32:40 +0200 Subject: [PATCH 14/17] Missing erase after remove_if --- src/Subsampling/benchmark/choose_n_farthest_points.cpp | 9 +++++++++ src/Subsampling/include/gudhi/choose_n_farthest_points.h | 3 ++- 2 files changed, 11 insertions(+), 1 deletion(-) 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 From 99eb8b7d26499ceef28e3b4e2032d8557d40e92f Mon Sep 17 00:00:00 2001 From: Marc Glisse Date: Wed, 20 Sep 2023 09:54:23 +0200 Subject: [PATCH 15/17] =?UTF-8?q?Cech=20->=20=C4=8Cech?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Vincent Rouvreau <10407034+VincentRouvreau@users.noreply.github.com> --- src/Cech_complex/include/gudhi/MEB_filtration.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Cech_complex/include/gudhi/MEB_filtration.h b/src/Cech_complex/include/gudhi/MEB_filtration.h index c1eefe8dbb..fd9837e873 100644 --- a/src/Cech_complex/include/gudhi/MEB_filtration.h +++ b/src/Cech_complex/include/gudhi/MEB_filtration.h @@ -21,7 +21,7 @@ namespace Gudhi::cech_complex { * each simplex a filtration value equal to the squared radius of its minimal * enclosing ball. * - * Applied on a Cech complex, it recomputes the same values (squared). Applied on a Delaunay triangulation, it computes the Delaunay-Cech filtration. + * 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`. From 3755f5d02289c2279ffaf231e428f8b7b4b91f08 Mon Sep 17 00:00:00 2001 From: Marc Glisse Date: Wed, 20 Sep 2023 09:54:53 +0200 Subject: [PATCH 16/17] MEB acronym Co-authored-by: Vincent Rouvreau <10407034+VincentRouvreau@users.noreply.github.com> --- src/Cech_complex/include/gudhi/MEB_filtration.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Cech_complex/include/gudhi/MEB_filtration.h b/src/Cech_complex/include/gudhi/MEB_filtration.h index fd9837e873..86d14df86d 100644 --- a/src/Cech_complex/include/gudhi/MEB_filtration.h +++ b/src/Cech_complex/include/gudhi/MEB_filtration.h @@ -19,7 +19,7 @@ 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 radius of its minimal - * enclosing ball. + * 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. * From 715c093450941aa82d2fe178f052c5321e542ae8 Mon Sep 17 00:00:00 2001 From: Vincent Rouvreau Date: Wed, 20 Sep 2023 11:08:19 +0200 Subject: [PATCH 17/17] warns if Hera submodule is not initialzied --- src/cmake/modules/GUDHI_submodules.cmake | 6 ++++++ 1 file changed, 6 insertions(+) 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