Merge branch 'master' into master

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_

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.

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<Point_d, FT>;
+  using Sphere = std::pair<Point_d, FT>;
 
  private:
 

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 <a href="https://doc.cgal.org/latest/Kernel_d/structCGAL_1_1Epeck__d.html">CGAL::Epeck_d</a>, the filtration values are computed exactly. Default is false.
+ */
+
+template<typename Kernel, typename SimplicialComplexForMEB, typename PointRange>
+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<Point_d, FT>;
+
+  using Vertex_handle = typename SimplicialComplexForMEB::Vertex_handle;
+  using Simplex_handle = typename SimplicialComplexForMEB::Simplex_handle;
+  using Filtration_value = typename SimplicialComplexForMEB::Filtration_value;
+
+  std::vector<Sphere> cache_;
+  std::vector<Point_d> pts;
+  CGAL::NT_converter<FT, Filtration_value> 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_

src/Cech_complex/test/test_cech_complex.cpp

@@ -19,6 +19,7 @@
 #include <algorithm>  // std::max
 
 #include <gudhi/Cech_complex.h>
+#include <gudhi/MEB_filtration.h>
 // to construct Cech_complex from a OFF file of points
 #include <gudhi/Points_off_io.h>
 #include <gudhi/Simplex_tree.h>
@@ -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) {

src/Collapse/include/gudhi/Flag_complex_edge_collapser.h

@@ -283,7 +283,7 @@ struct Flag_complex_edge_collapser {
 };
 
 template<class R> R to_range(R&& r) { return std::move(r); }
-template<class R, class T> R to_range(T&& t) { R r; r.insert(r.end(), t.begin(), t.end()); return r; }
+template<class R, class T> R to_range(T const& t) { R r; r.insert(r.end(), t.begin(), t.end()); return r; }
 
 template<class FilteredEdgeRange, class Delay>
 auto flag_complex_collapse_edges(FilteredEdgeRange&& edges, Delay&&delay) {

src/Subsampling/benchmark/choose_n_farthest_points.cpp

@@ -76,7 +76,16 @@ int main(int argc, char**argv) {
             << " vs metric "              << std::chrono::duration_cast<std::chrono::milliseconds>((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

src/Subsampling/include/gudhi/choose_n_farthest_points.h

@@ -320,14 +320,15 @@ 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.
             if (d <= 3 * radius) { l_neighbors.insert(near); }
             // 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

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()

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()