Euler characteristic + magnitude

src/Simplex_tree/include/gudhi/Simplex_tree.h

@@ -1701,6 +1701,71 @@ class Simplex_tree {
     }
   }
 
+ public:
+  /** \brief Returns a quantity similar to the Euler characteristic. It computes the image by `func`
+   * of the filtration value for each simplex, multiplies with \f$(-1)^{dim}\f$ and sums them all.
+   * The Euler characteristic corresponds to the case \f$func(f)=1\f$, while magnitude uses \f$func(f)=e^{-t\cdot f}\f$.
+   * @param[in] func function that takes a Filtration_value and returns a number.
+   */
+  template<class Func>
+  auto euler(Func&& func){
+    return rec_euler(&root_, func);
+  }
+  /** \brief Returns the Euler characteristic of the simplicial complex. */
+  size_t euler_characteristic(){
+    return euler([](auto){return 1;});
+  }
+  /** \brief Returns the magnitude of the filtration.
+   * @param[in] t parameter of the magnitude.
+   */
+  Filtration_value magnitude(Filtration_value t){
+    return euler([=](Filtration_value f){using std::exp; return exp(-t*f);});
+  }
+ private:
+  template<class F>
+  auto rec_euler(Siblings* sib, F& f) -> std::decay_t<decltype(f(std::declval<Filtration_value>()))> {
+    // We could use a style closer to visitors and pass the current value as argument of f.
+    std::decay_t<decltype(f(std::declval<Filtration_value>()))> ret = 0;
+    // TODO: use parallel_reduce or similar
+    for (auto sh = sib->members().begin(); sh != sib->members().end(); ++sh) {
+      ret += f(filtration(sh));
+      if (has_children(sh))
+        ret -= rec_euler(sh->second.children(), f);
+    }
+    return ret;
+  }
+
+// Print a Simplex_tree in os.
+  friend std::ostream& operator<<(std::ostream& os, Simplex_tree& st) {
+    for (auto sh : st.filtration_simplex_range()) {
+      os << st.dimension(sh) << " ";
+      for (auto v : st.simplex_vertex_range(sh)) {
+        os << v << " ";
+      }
+      os << st.filtration(sh) << "\n";  // TODO(VR): why adding the key ?? not read ?? << "     " << st.key(sh) << " \n";
+    }
+    return os;
+  }
+
+  friend std::istream& operator>>(std::istream& is, Simplex_tree& st) {
+    std::vector<Vertex_handle> simplex;
+    Filtration_value fil;
+    int max_dim = -1;
+    while (read_simplex(is, simplex, fil)) {
+      // read all simplices in the file as a list of vertices
+      // Warning : simplex_size needs to be casted in int - Can be 0
+      int dim = static_cast<int> (simplex.size() - 1);
+      if (max_dim < dim) {
+        max_dim = dim;
+      }
+      // insert every simplex in the simplex tree
+      st.insert_simplex(simplex, fil);
+      simplex.clear();
+    }
+    st.set_dimension(max_dim);
+    return is;
+  }
+
  private:
   Vertex_handle null_vertex_;
   /** \brief Total number of simplices in the complex, without the empty simplex.*/
@@ -1713,45 +1778,10 @@ class Simplex_tree {
   bool dimension_to_be_lowered_ = false;
 };
 
-// Print a Simplex_tree in os.
-template<typename...T>
-std::ostream& operator<<(std::ostream & os, Simplex_tree<T...> & st) {
-  for (auto sh : st.filtration_simplex_range()) {
-    os << st.dimension(sh) << " ";
-    for (auto v : st.simplex_vertex_range(sh)) {
-      os << v << " ";
-    }
-    os << st.filtration(sh) << "\n";  // TODO(VR): why adding the key ?? not read ?? << "     " << st.key(sh) << " \n";
-  }
-  return os;
-}
-
-template<typename...T>
-std::istream& operator>>(std::istream & is, Simplex_tree<T...> & st) {
-  typedef Simplex_tree<T...> ST;
-  std::vector<typename ST::Vertex_handle> simplex;
-  typename ST::Filtration_value fil;
-  int max_dim = -1;
-  while (read_simplex(is, simplex, fil)) {
-    // read all simplices in the file as a list of vertices
-    // Warning : simplex_size needs to be casted in int - Can be 0
-    int dim = static_cast<int> (simplex.size() - 1);
-    if (max_dim < dim) {
-      max_dim = dim;
-    }
-    // insert every simplex in the simplex tree
-    st.insert_simplex(simplex, fil);
-    simplex.clear();
-  }
-  st.set_dimension(max_dim);
-
-  return is;
-}
-
 /** Model of SimplexTreeOptions.
  * 
  * Maximum number of simplices to compute persistence is <CODE>std::numeric_limits<std::uint32_t>::max()</CODE>
- * (about 4 billions of simplices). */
+ * (about 4 billion simplices). */
 struct Simplex_tree_options_full_featured {
   typedef linear_indexing_tag Indexing_tag;
   typedef int Vertex_handle;
@@ -1766,7 +1796,7 @@ struct Simplex_tree_options_full_featured {
  * `contiguous_vertices` option.
  * 
  * Maximum number of simplices to compute persistence is <CODE>std::numeric_limits<std::uint32_t>::max()</CODE>
- * (about 4 billions of simplices). */
+ * (about 4 billion simplices). */
 
 struct Simplex_tree_options_fast_persistence {
   typedef linear_indexing_tag Indexing_tag;

src/Simplex_tree/test/CMakeLists.txt

@@ -40,3 +40,9 @@ if (TBB_FOUND)
   target_link_libraries(Simplex_tree_graph_expansion_test_unit ${TBB_LIBRARIES})
 endif()
 gudhi_add_boost_test(Simplex_tree_graph_expansion_test_unit)
+
+add_executable ( Simplex_tree_euler_test_unit simplex_tree_euler.cpp )
+if (TBB_FOUND)
+  target_link_libraries(Simplex_tree_euler_test_unit ${TBB_LIBRARIES})
+endif()
+gudhi_add_boost_test(Simplex_tree_euler_test_unit)

src/Simplex_tree/test/simplex_tree_euler.cpp

@@ -0,0 +1,41 @@
+/*    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) 2022 Inria
+ *
+ *    Modification(s):
+ *      - YYYY/MM Author: Description of the modification
+ */
+
+#include <iostream>
+#include <fstream>
+#include <string>
+#include <algorithm>
+#include <utility>  // std::pair, std::make_pair
+#include <cmath>  // float comparison
+#include <limits>
+#include <functional>  // greater
+#include <tuple>  // std::tie
+
+#define BOOST_TEST_DYN_LINK
+#define BOOST_TEST_MODULE "simplex_tree"
+#include <boost/test/unit_test.hpp>
+
+#include "gudhi/Simplex_tree.h"
+
+using namespace Gudhi;
+
+BOOST_AUTO_TEST_CASE(euler) {
+  typedef Simplex_tree<> ST;
+  ST st;
+  BOOST_CHECK(st.euler_characteristic() == 0);
+  BOOST_CHECK(st.magnitude(-.5) == 0.);
+  st.insert_simplex_and_subfaces({0,1,2},.5);
+  BOOST_CHECK(st.euler_characteristic() == 1);
+  BOOST_CHECK(std::abs(st.magnitude(0) - 1.) < .000001);
+  BOOST_CHECK(std::abs(st.magnitude(2.) - std::exp(-1.)) < .000001);
+  st.insert_simplex({0},0);
+  BOOST_CHECK(st.euler_characteristic() == 1);
+  BOOST_CHECK(std::abs(st.magnitude(2.) - 1.) < .000001);
+}

src/python/gudhi/simplex_tree.pxd

@@ -67,6 +67,9 @@ cdef extern from "Simplex_tree_interface.h" namespace "Gudhi":
         Simplex_tree_interface_full_featured* collapse_edges(int nb_collapse_iteration) nogil except +
         void reset_filtration(double filtration, int dimension) nogil
         bint operator==(Simplex_tree_interface_full_featured) nogil
+        double euler_wrap(double(*)(void*,double),void*)
+        int euler_characteristic() nogil
+        double magnitude(double) nogil
         # Iterators over Simplex tree
         pair[vector[int], double] get_simplex_and_filtration(Simplex_tree_simplex_handle f_simplex) nogil
         Simplex_tree_simplices_iterator get_simplices_iterator_begin() nogil

src/python/gudhi/simplex_tree.pyx

@@ -16,9 +16,12 @@ __author__ = "Vincent Rouvreau"
 __copyright__ = "Copyright (C) 2016 Inria"
 __license__ = "MIT"
 
-cdef bool callback(vector[int] simplex, void *blocker_func):
+cdef bool callback_blocker(vector[int] simplex, void *blocker_func):
     return (<object>blocker_func)(simplex)
 
+cdef double callback_euler(void* fun, double filt):
+    return (<object>fun)(filt)
+
 # SimplexTree python interface
 cdef class SimplexTree:
     """The simplex tree is an efficient and flexible data structure for
@@ -497,7 +500,7 @@ cdef class SimplexTree:
         :param blocker_func: Blocker oracle.
         :type blocker_func: Callable[[List[int]], bool]
         """
-        self.get_ptr().expansion_with_blockers_callback(max_dim, callback, <void*>blocker_func)
+        self.get_ptr().expansion_with_blockers_callback(max_dim, callback_blocker, <void*>blocker_func)
 
     def persistence(self, homology_coeff_field=11, min_persistence=0, persistence_dim_max = False):
         """This function computes and returns the persistence of the simplicial complex.
@@ -700,10 +703,47 @@ cdef class SimplexTree:
 
     def __eq__(self, other:SimplexTree):
         """Test for structural equality
+
         :returns: True if the 2 simplex trees are equal, False otherwise.
         :rtype: bool
         """
         return dereference(self.get_ptr()) == dereference(other.get_ptr())
 
+    def euler(self, fun):
+        """
+        Computes a quantity similar to the Euler characteristic. It computes the image by `fun`
+        of the filtration value for each simplex, negates it for odd dimensions and sums them all.
+        The Euler characteristic corresponds to the case `lambda _:1` while magnitude uses `lambda f:math.exp(-t*f)`.
+
+        :param fun: callback that computes some number from a filtration value.
+        :type fun: Callable[[float], float]
+        :returns: the quantity defined above.
+        :rtype: float
+        """
+        return self.get_ptr().euler_wrap(callback_euler, <void*>fun)
+
+    def euler_characteristic(self) -> int:
+        """
+        :returns: the Euler characteristic of the simplicial complex.
+        :rtype: int
+        """
+        cdef int ret
+        with nogil:
+            ret = self.get_ptr().euler_characteristic()
+        return ret
+
+    def magnitude(self, t: float) -> float:
+        """
+        :param t: parameter of the magnitude.
+        :type t: int
+        :returns: the magnitude of the filtration.
+        :rtype: float
+        """
+        cdef double t_ = t
+        cdef double ret
+        with nogil:
+            ret = self.get_ptr().magnitude(t_)
+        return ret
+
 cdef intptr_t _get_copy_intptr(SimplexTree stree) nogil:
     return <intptr_t>(new Simplex_tree_interface_full_featured(dereference(stree.get_ptr())))

src/python/include/Simplex_tree_interface.h

@@ -203,6 +203,10 @@ class Simplex_tree_interface : public Simplex_tree<SimplexTreeOptions> {
     });
   }
 
+  double euler_wrap(double (*call)(void*, double),void* fun) {
+    return Base::euler([=](double filt){ return call(fun, filt); });
+  }
+
   // Iterator over the simplex tree
   Complex_simplex_iterator get_simplices_iterator_begin() {
     // this specific case works because the range is just a pair of iterators - won't work if range was a vector

src/python/test/test_simplex_tree.py

@@ -11,6 +11,7 @@
 from gudhi import SimplexTree, __GUDHI_USE_EIGEN3
 import numpy as np
 import pytest
+import math
 
 __author__ = "Vincent Rouvreau"
 __copyright__ = "Copyright (C) 2016 Inria"
@@ -560,3 +561,16 @@ def blocker(simplex):
     assert st.filtration([0,2,3]) == 6.
     assert st.filtration([1,2,3]) == 6.
     assert st.filtration([0,1,2,3]) == 7.
+
+class Magn:
+    def __init__(self, tt):
+        self.t = tt
+    def __call__(self, f):
+        return math.exp(-f*self.t)
+
+def test_euler():
+    st = SimplexTree()
+    st.insert([0, 1], 2)
+    assert st.euler_characteristic() == 1
+    assert st.magnitude(-0.5) == pytest.approx(np.e)
+    assert st.euler(Magn(-0.5)) == pytest.approx(np.e)