From 36de6423be68cc113853eaf0388624d4b9666f7d Mon Sep 17 00:00:00 2001 From: David Loiseaux Date: Mon, 25 Sep 2023 15:28:42 +0200 Subject: [PATCH] safe+slow simplextree conversion and optimizations --- src/python/gudhi/simplex_tree_multi.pyx | 30 ++++++++-- .../include/Simplex_tree_interface_multi.h | 36 +----------- src/python/include/Simplex_tree_multi.h | 58 +++++++++++-------- 3 files changed, 60 insertions(+), 64 deletions(-) diff --git a/src/python/gudhi/simplex_tree_multi.pyx b/src/python/gudhi/simplex_tree_multi.pyx index c3f625033a..d2beda7069 100644 --- a/src/python/gudhi/simplex_tree_multi.pyx +++ b/src/python/gudhi/simplex_tree_multi.pyx @@ -59,8 +59,8 @@ cdef extern from "Simplex_tree_multi.h" namespace "Gudhi::multiparameter": vector[vector[vector[value_type]]] get_filtration_values(uintptr_t, const vector[int]&) except + nogil -cdef bool callback(vector[int] simplex, void *blocker_func): - return (blocker_func)(simplex) +# cdef bool callback(vector[int] simplex, void *blocker_func): + # return (blocker_func)(simplex) # SimplexTree python interface cdef class SimplexTreeMulti: @@ -84,7 +84,7 @@ cdef class SimplexTreeMulti: # cdef Simplex_tree_persistence_interface * pcohptr # Fake constructor that does nothing but documenting the constructor - def __init__(self, other = None, num_parameters:int=2,default_values=[]): + def __init__(self, other = None, num_parameters:int=2,default_values=[], safe_conversion=False): """SimplexTreeMulti constructor. :param other: If `other` is `None` (default value), an empty `SimplexTreeMulti` is created. @@ -102,6 +102,7 @@ cdef class SimplexTreeMulti: # The real cython constructor def __cinit__(self, other = None, int num_parameters=2, default_values=[-np.inf], # I'm not sure why `[]` does not work. Cython bug ? +bool safe_conversion=False, ): #TODO doc cdef vector[value_type] c_default_values=default_values if other is not None: @@ -111,7 +112,11 @@ cdef class SimplexTreeMulti: self.filtration_grid = other.filtration_grid elif isinstance(other, SimplexTree): # Constructs a SimplexTreeMulti from a SimplexTree self.thisptr = (new Simplex_tree_multi_interface()) - multify_from_ptr(other.thisptr, self.thisptr, num_parameters, c_default_values) + if safe_conversion: + new_st_multi = _safe_simplextree_multify(other, num_parameters = num_parameters) + self.thisptr, new_st_multi.thisptr = new_st_multi.thisptr, self.thisptr + else: + multify_from_ptr(other.thisptr, self.thisptr, num_parameters, c_default_values) else: raise TypeError("`other` argument requires to be of type `SimplexTree`, `SimplexTreeMulti`, or `None`.") else: @@ -1217,4 +1222,19 @@ def _simplextree_multify(simplextree:SimplexTree, num_parameters:int=2, default_ cdef vector[value_type] c_default_values=default_values with nogil: multify_from_ptr(old_ptr, new_ptr, c_num_parameters, c_default_values) - return \ No newline at end of file + return st + +def _safe_simplextree_multify(simplextree:SimplexTree,int num_parameters=2)->SimplexTreeMulti: + if isinstance(simplextree, SimplexTreeMulti): + return simplextree + simplices = [[] for _ in range(simplextree.dimension()+1)] + filtration_values = [[] for _ in range(simplextree.dimension()+1)] + for s,f in simplextree.get_simplices(): + filtration_values[len(s)-1].append(f) + simplices[len(s)-1].append(s) + st_multi = SimplexTreeMulti(num_parameters=1) + for batch_simplices, batch_filtrations in zip(simplices,filtration_values): + st_multi.insert_batch(np.asarray(batch_simplices).T, np.asarray(batch_filtrations)[:,None]) + if num_parameters > 1: + st_multi.set_num_parameter(num_parameters) + return st_multi \ No newline at end of file diff --git a/src/python/include/Simplex_tree_interface_multi.h b/src/python/include/Simplex_tree_interface_multi.h index db493f9302..139def9ee2 100644 --- a/src/python/include/Simplex_tree_interface_multi.h +++ b/src/python/include/Simplex_tree_interface_multi.h @@ -212,40 +212,6 @@ class Simplex_tree_interface_multi : public Simplex_tree_interface> &point_list){ // TODO multi-critical -// const int npts = point_list.size(); -// if (npts == 0){ -// return {}; -// } -// using Gudhi::multi_filtrations::Finitely_critical_multi_filtration; - -// euler_chars_type out(point_list.size(), 0.); - -// // auto is_greater = [nparameters](const point_type &a, const point_type &b){ //french greater -// // for (int i = 0; i< nparameters; i++) -// // if( a[i] < b[i]) -// // return false; -// // return true; -// // }; -// // #pragma omp parallel for -// for (int i = 0; i< npts; i++){ // Maybe add a pragma here for parallel -// auto &euler_char_at_point = out[i]; -// // #pragma omp parallel for reduction(+:euler_char_at_point) // GUDHI : not possible, need a RANDOM ACCESS ITERATOR -// for(const auto &SimplexHandle : Base::complex_simplex_range()){ -// // const Finitely_critical_multi_filtration &pt = *(Finitely_critical_multi_filtration*)(&point_list[i]); -// options_multi::Filtration_value filtration = Base::filtration(SimplexHandle); -// // const Finitely_critical_multi_filtration &filtration = *(Finitely_critical_multi_filtration*)(&filtration_); -// Finitely_critical_multi_filtration pt(point_list[i]); -// // if (is_greater(pt, filtration)){ -// if (filtration <= pt){ -// int sign = Base::dimension(SimplexHandle) %2 ? -1 : 1; -// euler_char_at_point += sign; -// } -// } -// } -// return out; -// } void resize_all_filtrations(int num){ //TODO : that is for 1 critical filtrations if (num < 0) return; for(const auto &SimplexHandle : Base::complex_simplex_range()){ @@ -260,7 +226,7 @@ class Simplex_tree_interface_multi : public Simplex_tree_interface; +using interface_std = Simplex_tree; // Interface not necessary (smaller so should do less segfaults) using interface_multi = Simplex_tree_interface_multi; diff --git a/src/python/include/Simplex_tree_multi.h b/src/python/include/Simplex_tree_multi.h index 941298ee53..2fcdf35821 100644 --- a/src/python/include/Simplex_tree_multi.h +++ b/src/python/include/Simplex_tree_multi.h @@ -61,9 +61,7 @@ simplextreeinterface& get_simplextree_from_pointer(const uintptr_t splxptr){ //D } template void multify(simplextree_std &st, simplextree_multi &st_multi, const int num_parameters, const typename simplextree_multi::Options::Filtration_value& default_values={}){ - if (num_parameters <= 0) - {std::cerr << "Empty filtration\n"; throw ;} - // if (default_values.size() -1 > num_parameters) + // if (default_values.size() -1 > num_parameters) // {std::cerr << "default values too large !\n"; throw ;} typename simplextree_multi::Options::Filtration_value f(num_parameters); for (auto i = 0u; i(default_values.size()), static_cast(num_parameters-1));i++) @@ -74,9 +72,12 @@ void multify(simplextree_std &st, simplextree_multi &st_multi, const int num_par simplex.clear(); for (auto vertex : st.simplex_vertex_range(simplex_handle)) simplex.push_back(vertex); + + if (num_parameters > 0) f[0] = st.filtration(simplex_handle); // std::cout << "Inserting " << f << "\n"; - st_multi.insert_simplex(simplex,f); + auto filtration_copy = f; + st_multi.insert_simplex(simplex,filtration_copy); } } @@ -130,18 +131,28 @@ void flatten_diag(simplextree_std &st, simplextree_multi &st_multi, const std::v -/// @brief projects the point x on the grid +/// @brief turns filtration value x into coordinates in the grid /// @tparam out_type /// @param x /// @param grid /// @return -template -std::vector find_coordinates(const std::vector &x, const multi_filtration_grid &grid){ +template +inline void find_coordinates(vector_like& x, const multi_filtration_grid &grid){ // TODO: optimize with, e.g., dichotomy - std::vector coordinates(grid.size()); - for (int parameter = 0; parameter < (int)grid.size(); parameter++){ + + for (auto parameter = 0u; parameter < grid.size(); parameter++){ const auto& filtration = grid[parameter]; // assumes its sorted const auto to_project = x[parameter]; + if constexpr (std::numeric_limits::has_infinity) + if (to_project == std::numeric_limits::infinity()){ + x[parameter] = std::numeric_limits::infinity(); +continue; + } + + if (to_project >= filtration.back()){ + x[parameter] = filtration.size()-1; + continue; + } // deals with infinite value at the end of the grid, TODO DEAL unsigned int i = 0; // std::cout << to_project<< " " < filtration[i] && i find_coordinates(const std::vector distance_vector(filtration.size()); // for (int i = 0; i < (int)filtration.size(); i++){ @@ -162,29 +173,26 @@ std::vector find_coordinates(const std::vector &st_multi = *(Gudhi::Simplex_tree*)(splxptr); - int num_parameters = st_multi.get_number_of_parameters(); - if ((int)grid.size() != num_parameters){ + auto num_parameters = static_cast(st_multi.get_number_of_parameters()); + if (grid.size() != num_parameters){ std::cerr << "Bad grid !" << std::endl; throw; } for (const auto &simplex_handle : st_multi.complex_simplex_range()){ - std::vector simplex_filtration = st_multi.filtration(simplex_handle); - if (coordinate_values) - st_multi.assign_filtration(simplex_handle, find_coordinates(simplex_filtration, grid)); - else{ - auto coordinates = find_coordinates(simplex_filtration, grid); - std::vector squeezed_filtration(num_parameters); - for (int parameter = 0; parameter < num_parameters; parameter++) - squeezed_filtration[parameter] = grid[parameter][coordinates[parameter]]; - st_multi.assign_filtration(simplex_handle, squeezed_filtration); + auto& simplex_filtration = st_multi.filtration_mutable(simplex_handle); + find_coordinates(simplex_filtration, grid); // turns the simplexfiltration into coords in the grid + if (!coordinate_values){ + for (auto parameter = 0u; parameter < num_parameters; parameter++) + simplex_filtration[parameter] = grid[parameter][simplex_filtration[parameter]]; } } return; @@ -214,6 +222,8 @@ std::vector get_filtration_values(const uintptr_t splxptr } // namespace Gudhi + + namespace std { template<>