From 20e53d0016391aaa3d980d176aabb426043f94ce Mon Sep 17 00:00:00 2001 From: Marc Glisse Date: Sat, 14 May 2022 23:49:17 +0200 Subject: [PATCH 1/3] Initial implementation of Euler characteristic --- src/Simplex_tree/include/gudhi/Simplex_tree.h | 104 +++++++++++------- src/Simplex_tree/test/CMakeLists.txt | 6 + src/Simplex_tree/test/simplex_tree_euler.cpp | 41 +++++++ 3 files changed, 114 insertions(+), 37 deletions(-) create mode 100644 src/Simplex_tree/test/simplex_tree_euler.cpp diff --git a/src/Simplex_tree/include/gudhi/Simplex_tree.h b/src/Simplex_tree/include/gudhi/Simplex_tree.h index d48f6616d5..d251374168 100644 --- a/src/Simplex_tree/include/gudhi/Simplex_tree.h +++ b/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 + 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 + auto rec_euler(Siblings* sib, F& f) -> std::decay_t()))> { + // We could use a style closer to visitors and pass the current value as argument of f. + std::decay_t()))> 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 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 (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 -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; -} - -template -std::istream& operator>>(std::istream & is, Simplex_tree & st) { - typedef Simplex_tree ST; - std::vector 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 (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 std::numeric_limits::max() - * (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 std::numeric_limits::max() - * (about 4 billions of simplices). */ + * (about 4 billion simplices). */ struct Simplex_tree_options_fast_persistence { typedef linear_indexing_tag Indexing_tag; diff --git a/src/Simplex_tree/test/CMakeLists.txt b/src/Simplex_tree/test/CMakeLists.txt index 25b562e0db..ea50523daa 100644 --- a/src/Simplex_tree/test/CMakeLists.txt +++ b/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) diff --git a/src/Simplex_tree/test/simplex_tree_euler.cpp b/src/Simplex_tree/test/simplex_tree_euler.cpp new file mode 100644 index 0000000000..6b51bd20a0 --- /dev/null +++ b/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 +#include +#include +#include +#include // std::pair, std::make_pair +#include // float comparison +#include +#include // greater +#include // std::tie + +#define BOOST_TEST_DYN_LINK +#define BOOST_TEST_MODULE "simplex_tree" +#include + +#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); +} From 97386c209d1837cecc946fa17cbda48d8aff161f Mon Sep 17 00:00:00 2001 From: Marc Glisse Date: Sun, 15 May 2022 00:28:45 +0200 Subject: [PATCH 2/3] Python interface to Euler, but not the generic version --- src/python/gudhi/simplex_tree.pxd | 2 ++ src/python/gudhi/simplex_tree.pyx | 23 +++++++++++++++++++++++ src/python/test/test_simplex_tree.py | 6 ++++++ 3 files changed, 31 insertions(+) diff --git a/src/python/gudhi/simplex_tree.pxd b/src/python/gudhi/simplex_tree.pxd index 4f22966342..aa34f29493 100644 --- a/src/python/gudhi/simplex_tree.pxd +++ b/src/python/gudhi/simplex_tree.pxd @@ -67,6 +67,8 @@ 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 + 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 diff --git a/src/python/gudhi/simplex_tree.pyx b/src/python/gudhi/simplex_tree.pyx index a491418443..d726745ffc 100644 --- a/src/python/gudhi/simplex_tree.pyx +++ b/src/python/gudhi/simplex_tree.pyx @@ -705,5 +705,28 @@ cdef class SimplexTree: """ return dereference(self.get_ptr()) == dereference(other.get_ptr()) + 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 (new Simplex_tree_interface_full_featured(dereference(stree.get_ptr()))) diff --git a/src/python/test/test_simplex_tree.py b/src/python/test/test_simplex_tree.py index eb481a49b0..bf5102cae3 100755 --- a/src/python/test/test_simplex_tree.py +++ b/src/python/test/test_simplex_tree.py @@ -560,3 +560,9 @@ 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. + +def test_euler(): + st = SimplexTree() + st.insert([0, 1], 2) + assert st.euler_characteristic() == 1 + assert st.magnitude(-0.5) == pytest.approx(np.e) From 0d2f94f5b68475a5512522c610dc865c152d4352 Mon Sep 17 00:00:00 2001 From: Marc Glisse Date: Sun, 15 May 2022 17:23:44 +0200 Subject: [PATCH 3/3] Python wrapper for generic Euler-type functions --- src/python/gudhi/simplex_tree.pxd | 1 + src/python/gudhi/simplex_tree.pyx | 21 +++++++++++++++++++-- src/python/include/Simplex_tree_interface.h | 4 ++++ src/python/test/test_simplex_tree.py | 8 ++++++++ 4 files changed, 32 insertions(+), 2 deletions(-) diff --git a/src/python/gudhi/simplex_tree.pxd b/src/python/gudhi/simplex_tree.pxd index aa34f29493..94b4980101 100644 --- a/src/python/gudhi/simplex_tree.pxd +++ b/src/python/gudhi/simplex_tree.pxd @@ -67,6 +67,7 @@ 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 diff --git a/src/python/gudhi/simplex_tree.pyx b/src/python/gudhi/simplex_tree.pyx index d726745ffc..b97f1e8997 100644 --- a/src/python/gudhi/simplex_tree.pyx +++ b/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 (blocker_func)(simplex) +cdef double callback_euler(void* fun, double filt): + return (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, blocker_func) + self.get_ptr().expansion_with_blockers_callback(max_dim, callback_blocker, 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,11 +703,25 @@ 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, fun) + def euler_characteristic(self) -> int: """ :returns: the Euler characteristic of the simplicial complex. diff --git a/src/python/include/Simplex_tree_interface.h b/src/python/include/Simplex_tree_interface.h index aa3dac189b..af627d591d 100644 --- a/src/python/include/Simplex_tree_interface.h +++ b/src/python/include/Simplex_tree_interface.h @@ -203,6 +203,10 @@ class Simplex_tree_interface : public Simplex_tree { }); } + 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 diff --git a/src/python/test/test_simplex_tree.py b/src/python/test/test_simplex_tree.py index bf5102cae3..b7e518d279 100755 --- a/src/python/test/test_simplex_tree.py +++ b/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" @@ -561,8 +562,15 @@ def blocker(simplex): 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)