diff --git a/src/Simplex_tree/concept/FiltrationValue.h b/src/Simplex_tree/concept/FiltrationValue.h index dd06b5d8e6..9c71778ab2 100644 --- a/src/Simplex_tree/concept/FiltrationValue.h +++ b/src/Simplex_tree/concept/FiltrationValue.h @@ -11,9 +11,12 @@ /** \brief Value type for a filtration function on a cell complex. * - * Needs to implement `std::numeric_limits::has_infinity` and when it returns `true`, - * has to implement `std::numeric_limits::infinity()`. If it returns `false`, has to - * implement `std::numeric_limits::max()` instead. + * Needs to implement `std::numeric_limits::has_infinity`, + * `std::numeric_limits::infinity()` and `std::numeric_limits::max()`. + * But when `std::numeric_limits::has_infinity` returns `true`, + * `std::numeric_limits::max()` can simply throw when called, as well as, + * `std::numeric_limits::infinity()` if `std::numeric_limits::has_infinity` + * returns `false`. * * A filtration of a cell complex (see FilteredComplex) is * a function \f$f:\mathbf{K} \rightarrow \mathbb{R}\f$ satisfying \f$f(\tau)\leq @@ -32,38 +35,42 @@ struct FiltrationValue { /** * @brief Strictly smaller operator. If the filtration values are totally ordered, should be a StrictWeakOrdering. */ - bool operator<(FiltrationValue f1, FiltrationValue f2); + friend bool operator<(const FiltrationValue& f1, const FiltrationValue& f2); // only for prune_above_filtration /** * @brief Smaller or equal operator. */ - bool operator<=(FiltrationValue f1, FiltrationValue f2); + friend bool operator<=(const FiltrationValue& f1, const FiltrationValue& f2); /** * @brief Equality operator */ - bool operator==(FiltrationValue f1, FiltrationValue f2); + friend bool operator==(const FiltrationValue& f1, const FiltrationValue& f2); /** * @brief Not equal operator */ - bool operator!=(FiltrationValue f1, FiltrationValue f2); + friend bool operator!=(const FiltrationValue& f1, const FiltrationValue& f2); /** * @brief Given two filtration values at which a simplex exists, returns the minimal union of births generating * a lifetime including those two values. The overload for native arithmetic types like `double` or `int` is * already implemented. - * If the filtration is 1-critical and totally ordered (as in one parameter persistence), the union is simply - * the minimum of the two values. If the filtration is 1-critical, but not totally ordered (possible for - * multi-persistence), than the union is also the minium if the two given values are comparable and the - * method should throw an error if they are not, as a same simplex should not exist at those two values. - * Finally, if the filtration is k-critical, FiltrationValue should be able to store an union of values and this - * method adds the values of @p f2 in @p f1 and removes the values from @p f1 which are comparable and greater than - * other values. + * + * For a k-critical filtration, FiltrationValue should be able to store an union of values (corresponding to the + * different births of a same simplex) and this method adds the values of @p f2 in @p f1 and removes the values + * from @p f1 which are comparable and greater than other values. + * In the special case of 1-critical filtration, as the union should not contain more than one birth element, + * this method is expected to throw if the two given elements in the filtration values are not comparable. + * If they are comparable, the union is simply the minimum of both. + * + * @return True if and only if the values in @p f1 were actually modified. */ - friend void unify_births(FiltrationValue& f1, const FiltrationValue& f2); + friend bool unify_births(FiltrationValue& f1, const FiltrationValue& f2); /** * @brief Given two filtration values, stores in the first value the greatest common upper bound of the two values. * The overload for native arithmetic types like `double` or `int` is already implemented. + * + * @return True if and only if the values in @p f1 were actually modified. */ - friend void push_to_greatest_common_upper_bound(FiltrationValue& f1, const FiltrationValue& f2); + friend bool push_to_smallest_common_upper_bound(FiltrationValue& f1, const FiltrationValue& f2); }; diff --git a/src/Simplex_tree/include/gudhi/Simplex_tree.h b/src/Simplex_tree/include/gudhi/Simplex_tree.h index fefedd62ff..bde48e1e04 100644 --- a/src/Simplex_tree/include/gudhi/Simplex_tree.h +++ b/src/Simplex_tree/include/gudhi/Simplex_tree.h @@ -900,14 +900,12 @@ class Simplex_tree { if (!res_insert.second) { if constexpr (Options::store_filtration){ // if already in the complex - auto oldFil = res_insert.first->second.filtration(); - unify_births(res_insert.first->second.filtration_raw(), filtration); - if (oldFil == res_insert.first->second.filtration()) { - // if filtration value unchanged - return std::pair(null_simplex(), false); + if (unify_births(res_insert.first->second.filtration_raw(), filtration)) { + // if filtration value modified + return res_insert; } - // if filtration value modified - return res_insert; + // if filtration value unchanged + return std::pair(null_simplex(), false); } else { return std::pair(null_simplex(), false); } @@ -1030,9 +1028,7 @@ class Simplex_tree { update_simplex_tree_after_node_insertion(insertion_result.first); } else { if constexpr (Options::store_filtration){ - auto oldFil = simplex_one->second.filtration(); - unify_births(simplex_one->second.filtration_raw(), filt); - if (oldFil == simplex_one->second.filtration()) { + if (!unify_births(simplex_one->second.filtration_raw(), filt)) { // FIXME: this interface makes no sense, and it doesn't seem to be tested. insertion_result.first = null_simplex(); } @@ -1789,8 +1785,8 @@ class Simplex_tree { intersection.emplace_back(begin1->first, Node(nullptr, filtration_)); } else { Filtration_value filt = begin1->second.filtration(); - push_to_greatest_common_upper_bound(filt, begin2->second.filtration()); - push_to_greatest_common_upper_bound(filt, filtration_); + push_to_smallest_common_upper_bound(filt, begin2->second.filtration()); + push_to_smallest_common_upper_bound(filt, filtration_); intersection.emplace_back(begin1->first, Node(nullptr, filt)); } if (++begin1 == end1 || ++begin2 == end2) @@ -1859,7 +1855,7 @@ class Simplex_tree { to_be_inserted=false; break; } - push_to_greatest_common_upper_bound(filt, filtration(border_child)); + push_to_smallest_common_upper_bound(filt, filtration(border_child)); } if (to_be_inserted) { intersection.emplace_back(next->first, Node(nullptr, filt)); @@ -1993,18 +1989,13 @@ class Simplex_tree { auto fun = [&modified, this](Simplex_handle sh, int dim) -> void { if (dim == 0) return; - const Filtration_value_rep old = sh->second.filtration(); Filtration_value& max_filt_border_value = sh->second.filtration_raw(); // Find the maximum filtration value in the border and assigns it if it is greater than the current for (Simplex_handle b : boundary_simplex_range(sh)) { // considers NaN as the lowest possible value. - push_to_greatest_common_upper_bound(max_filt_border_value, b->second.filtration()); - } - - if (old != sh->second.filtration()) { - // Store the filtration modification information - modified = true; + // stores the filtration modification information + modified = push_to_smallest_common_upper_bound(max_filt_border_value, b->second.filtration()) || modified; } }; // Loop must be from the end to the beginning, as higher dimension simplex are always on the left part of the tree diff --git a/src/Simplex_tree/include/gudhi/Simplex_tree/simplex_tree_options.h b/src/Simplex_tree/include/gudhi/Simplex_tree/simplex_tree_options.h index 804f1bfc4b..41e528c37b 100644 --- a/src/Simplex_tree/include/gudhi/Simplex_tree/simplex_tree_options.h +++ b/src/Simplex_tree/include/gudhi/Simplex_tree/simplex_tree_options.h @@ -108,48 +108,47 @@ struct Get_simplex_data_type> { typedef * This is the overload for when `Filtration_value` is a native arithmetic type, like double, int etc. * Because the filtration values are totally ordered then, the union is simply the minimum of the two values. */ -template -std::enable_if_t> unify_births(Arithmetic_filtration_value& f1, - Arithmetic_filtration_value f2) +template > > +bool unify_births(Arithmetic_filtration_value& f1, Arithmetic_filtration_value f2) { - f1 = std::min(f1, f2); -} - -/** - * @private - * @brief Given two filtration values, stores in the first value the greatest common upper bound of the two values. - * If a filtration value has value `NaN`, it should be considered as the lowest value possible. - * This is the overload for when `Filtration_value` is a native floating type, like double, float etc. - * Because the filtration values are totally ordered then, the upper bound is always the maximum of the two values. - */ -template -std::enable_if_t > push_to_greatest_common_upper_bound( - Floating_filtration_value& f1, - Floating_filtration_value f2) -{ - if (std::isnan(f1)) { + if (f1 > f2){ f1 = f2; - return; + return true; } - - // Computes the max while handling NaN as lowest value. - f1 = !(f1 <= f2) ? f1 : f2; + return false; } /** * @private * @brief Given two filtration values, stores in the first value the greatest common upper bound of the two values. * If a filtration value has value `NaN`, it should be considered as the lowest value possible. - * This is the overload for when `Filtration_value` is a native integer type, like int, unsigned int etc. + * This is the overload for when `Filtration_value` is a native arithmetic type, like double, float, int etc. * Because the filtration values are totally ordered then, the upper bound is always the maximum of the two values. */ -template -std::enable_if_t > push_to_greatest_common_upper_bound( - Integral_filtration_value& f1, - Integral_filtration_value f2) +template > > +bool push_to_smallest_common_upper_bound(Floating_filtration_value& f1, Floating_filtration_value f2) { - // NaN not possible. - f1 = f1 < f2 ? f2 : f1; + if constexpr (std::is_floating_point_v) { + if (std::isnan(f1)) { + f1 = f2; + return !std::isnan(f1); + } + + // Computes the max while handling NaN as lowest value. + if (f1 == f2 || !(f1 <= f2)) return false; + + f1 = f2; + return true; + } else { + // NaN not possible. + if (f1 < f2){ + f1 = f2; + return true; + } + return false; + } } } // namespace Gudhi