Skip to content

Commit

Permalink
safe+slow simplextree conversion and optimizations
Browse files Browse the repository at this point in the history
  • Loading branch information
@DavidLapous
DavidLapous committed Sep 25, 2023
1 parent 32a6a49 commit 36de642
Show file tree
Hide file tree
Showing 3 changed files with 60 additions and 64 deletions.
30 changes: 25 additions & 5 deletions src/python/gudhi/simplex_tree_multi.pyx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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 (<object>blocker_func)(simplex)
# cdef bool callback(vector[int] simplex, void *blocker_func):
# return (<object>blocker_func)(simplex)

# SimplexTree python interface
cdef class SimplexTreeMulti:
Expand All @@ -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.
Expand All @@ -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:
Expand All @@ -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 = <intptr_t>(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:
Expand Down Expand Up @@ -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
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
36 changes: 1 addition & 35 deletions src/python/include/Simplex_tree_interface_multi.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -212,40 +212,6 @@ class Simplex_tree_interface_multi : public Simplex_tree_interface<Simplex_tree_
/* simplex_list.shrink_to_fit();*/
return simplex_list;
}
// DEPRECATED, USE COORDINATE SIMPLEX TREE
// euler_chars_type euler_char(const std::vector<std::vector<options_multi::value_type>> &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<options_multi::value_type> &pt = *(Finitely_critical_multi_filtration<options_multi::value_type>*)(&point_list[i]);
// options_multi::Filtration_value filtration = Base::filtration(SimplexHandle);
// // const Finitely_critical_multi_filtration<options_multi::value_type> &filtration = *(Finitely_critical_multi_filtration<options_multi::value_type>*)(&filtration_);
// Finitely_critical_multi_filtration<options_multi::value_type> 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()){
Expand All @@ -260,7 +226,7 @@ class Simplex_tree_interface_multi : public Simplex_tree_interface<Simplex_tree_
};


using interface_std = Simplex_tree_interface<Simplex_tree_options_full_featured>;
using interface_std = Simplex_tree<Simplex_tree_options_full_featured>; // Interface not necessary (smaller so should do less segfaults)
using interface_multi = Simplex_tree_interface_multi<Simplex_tree_options_multidimensional_filtration>;


Expand Down
58 changes: 34 additions & 24 deletions src/python/include/Simplex_tree_multi.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -61,9 +61,7 @@ simplextreeinterface& get_simplextree_from_pointer(const uintptr_t splxptr){ //D
}
template<class simplextree_std, class simplextree_multi>
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<std::min(static_cast<unsigned int>(default_values.size()), static_cast<unsigned int>(num_parameters-1));i++)
Expand All @@ -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);

}
}
Expand Down Expand Up @@ -130,61 +131,68 @@ 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<typename out_type=int>
std::vector<out_type> find_coordinates(const std::vector<options_multi::value_type> &x, const multi_filtration_grid &grid){
template<typename out_type=int, typename vector_like>
inline void find_coordinates(vector_like& x, const multi_filtration_grid &grid){
// TODO: optimize with, e.g., dichotomy
std::vector<out_type> 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<typename vector_like::value_type>::has_infinity)
if (to_project == std::numeric_limits<typename vector_like::value_type>::infinity()){
x[parameter] = std::numeric_limits<typename vector_like::value_type>::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 << "\n";
while (to_project > filtration[i] && i<filtration.size()) {
i++;
// std::cout << to_project<< " " <<filtration[i]<<" " <<i << "\n";
}
if (i==0)
coordinates[parameter] = 0;
x[parameter] = 0;
else if (i < filtration.size()){
options_multi::value_type d1,d2;
typename vector_like::value_type d1,d2;
d1 = std::abs(filtration[i-1] - to_project);
d2 = std::abs(filtration[i] - to_project);
coordinates[parameter] = d1<d2 ? i-1 : i;
x[parameter] = d1<d2 ? i-1 : i;
}
// std::vector<options_multi::value_type> distance_vector(filtration.size());
// for (int i = 0; i < (int)filtration.size(); i++){
// distance_vector[i] = std::abs(to_project - filtration[i]);
// }
// coordinates[parameter] = std::distance(distance_vector.begin(), std::min_element(distance_vector.begin(), distance_vector.end()));
}
return coordinates;
// return x;
}


// TODO integer filtrations, does this help with performance ?
// projects filtrations values to the grid. If coordinate_values is set to true, the filtration values are the coordinates of this grid
void squeeze_filtration(uintptr_t splxptr, const multi_filtration_grid &grid, bool coordinate_values=true){

Simplex_tree<options_multi> &st_multi = *(Gudhi::Simplex_tree<options_multi>*)(splxptr);
int num_parameters = st_multi.get_number_of_parameters();
if ((int)grid.size() != num_parameters){
auto num_parameters = static_cast<unsigned int>(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<options_multi::value_type> simplex_filtration = st_multi.filtration(simplex_handle);
if (coordinate_values)
st_multi.assign_filtration(simplex_handle, find_coordinates<options_multi::value_type>(simplex_filtration, grid));
else{
auto coordinates = find_coordinates<int>(simplex_filtration, grid);
std::vector<options_multi::value_type> 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<options_multi::value_type>(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;
Expand Down Expand Up @@ -214,6 +222,8 @@ std::vector<multi_filtration_grid> get_filtration_values(const uintptr_t splxptr

} // namespace Gudhi



namespace std {

template<>
Expand Down

0 comments on commit 36de642

Please sign in to comment.