diff --git a/CMakeLists.txt b/CMakeLists.txt index a0752b637a..714cfe484f 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -48,6 +48,8 @@ add_gudhi_module(Witness_complex) add_gudhi_module(Nerve_GIC) add_gudhi_module(Persistence_matrix) add_gudhi_module(Zigzag_persistence) +add_gudhi_module(Multi_persistence) +add_gudhi_module(Multi_filtration) # Include module CMake subdirectories # GUDHI_SUB_DIRECTORIES is managed in CMAKE_MODULE_PATH/GUDHI_modules.cmake diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index bed9459fc2..33755176bb 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -47,6 +47,8 @@ add_gudhi_module(Witness_complex) add_gudhi_module(Nerve_GIC) add_gudhi_module(Persistence_matrix) add_gudhi_module(Zigzag_persistence) +add_gudhi_module(Multi_persistence) +add_gudhi_module(Multi_filtration) set(GUDHI_BIBLIO_DIR ${CMAKE_SOURCE_DIR}) # For "make doxygen" - Requires GUDHI_USER_VERSION_DIR to be set diff --git a/src/Multi_filtration/doc/COPYRIGHT b/src/Multi_filtration/doc/COPYRIGHT new file mode 100644 index 0000000000..70cac8b762 --- /dev/null +++ b/src/Multi_filtration/doc/COPYRIGHT @@ -0,0 +1,12 @@ +The files of this directory are 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): Hannah Schreiber, David Loiseaux + +Copyright (C) 2024 Inria + +This gives everyone the freedoms to use openFrameworks in any context: +commercial or non-commercial, public or private, open or closed source. + +You should have received a copy of the MIT License along with this program. +If not, see https://opensource.org/licenses/MIT. diff --git a/src/Multi_filtration/doc/Intro_multi_filtration_values.h b/src/Multi_filtration/doc/Intro_multi_filtration_values.h new file mode 100644 index 0000000000..b2c1d6e562 --- /dev/null +++ b/src/Multi_filtration/doc/Intro_multi_filtration_values.h @@ -0,0 +1,36 @@ +/* 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): Hannah Schreiber + * + * Copyright (C) 2024 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#ifndef DOC_MULTI_FILTRATION_INTRO_H_ +#define DOC_MULTI_FILTRATION_INTRO_H_ + +// needs namespace for Doxygen to link on classes +namespace Gudhi { +namespace multi_filtration { + +/** \defgroup multi_filtration Multi-parameter Filtration Values + * @{ + * \author David Loiseaux + * + * \section multifiltrationintro Multi-parameter Filtration Values + * + * TODO: introduction to @ref One_critical_filtration and @ref Multi_critical_filtration + * + * \subsection multifiltrationexamples Examples + * + * Here is a list of examples using the module: + * \li TODO?. + * + * @} + */ +} // namespace persistence_fields +} // namespace Gudhi + +#endif // DOC_MULTI_FILTRATION_INTRO_H_ diff --git a/src/Multi_filtration/doc/radius_density_complex_without_interval.png b/src/Multi_filtration/doc/radius_density_complex_without_interval.png new file mode 100644 index 0000000000..c371900441 Binary files /dev/null and b/src/Multi_filtration/doc/radius_density_complex_without_interval.png differ diff --git a/src/Multi_filtration/include/gudhi/Multi_critical_filtration.h b/src/Multi_filtration/include/gudhi/Multi_critical_filtration.h new file mode 100644 index 0000000000..5946690962 --- /dev/null +++ b/src/Multi_filtration/include/gudhi/Multi_critical_filtration.h @@ -0,0 +1,951 @@ +/* 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): David Loiseaux + * + * Copyright (C) 2023 Inria + * + * Modification(s): + * - 2024/08 Hannah Schreiber: Optimization and correction + numeric_limits + doc + * - YYYY/MM Author: Description of the modification + */ + +/** + * @file Multi_critical_filtration.h + * @author David Loiseaux + * @brief Contains the @ref Gudhi::multi_filtration::Multi_critical_filtration class. + */ + +#ifndef MULTI_CRITICAL_FILTRATIONS_H_ +#define MULTI_CRITICAL_FILTRATIONS_H_ + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +namespace Gudhi::multi_filtration { + +/** + * @class Multi_critical_filtration multi_critical_filtration.h gudhi/multi_critical_filtration.h + * @ingroup multi_filtration + * + * @brief Class encoding the different generators, i.e., apparition times, of a \f$k\f$-critical + * \f$\mathbb R^n\f$-filtration value, e.g., the filtration of a simplex, or of the algebraic generator of a module + * presentation. The class can be used as a vector whose indices correspond each to a generator, i.e., a one-critical + * filtration value. Then, the indices of each generator correspond to a particular parameter. + * E.g., \f$ f[i][p] \f$ will be \f$ p^{\textit{th}} \f$ parameter of the \f$ i^{\textit{th}} \f$ generator + * of this filtration value with @ref Multi_critical_filtration \f$ f \f$. + * + * @details Overloads `std::numeric_limits` such that: + * - `std::numeric_limits >::has_infinity` returns `true`, + * - `std::numeric_limits >::infinity()` returns + * @ref Multi_critical_filtration::inf() "", + * - `std::numeric_limits >::minus_infinity()` returns + * @ref Multi_critical_filtration::minus_inf() "", + * - `std::numeric_limits >::max()` throws, + * - `std::numeric_limits >::max(g,n)` returns a @ref Multi_critical_filtration + * with `g` generators of `n` parameters evaluated at value `std::numeric_limits::max()`, + * - `std::numeric_limits >::quiet_NaN()` returns + * @ref Multi_critical_filtration::nan() "". + * + * Multi-critical filtrations are filtrations such that the lifetime of each object is union of positive cones in + * \f$\mathbb R^n\f$, e.g., + * - \f$ \{ x \in \mathbb R^2 : x \ge (1,2)\} \cap \{ x \in \mathbb R^2 : x \ge (2,1)\} \f$ is finitely critical, + * and more particularly 2-critical, while + * - \f$ \{ x \in \mathbb R^2 : x \ge \mathrm{epigraph}(y \mapsto e^{-y})\} \f$ is not. + * + * The particular case of 1-critical filtrations is handled by @ref One_critical_filtration "". + * + * @tparam T Arithmetic type of an entry for one parameter of a filtration value. Has to be **signed** and + * to implement `std::isnan(T)`, `std::numeric_limits::has_quiet_NaN`, `std::numeric_limits::quiet_NaN()`, + * `std::numeric_limits::has_infinity`, `std::numeric_limits::infinity()` and `std::numeric_limits::max()`. + * If `std::numeric_limits::has_infinity` returns `false`, a call to `std::numeric_limits::infinity()` + * can simply throw. Examples are the native types `double`, `float` and `int`. + * @tparam co If `true`, reverses the poset order, i.e., the order \f$ \le \f$ in \f$ \mathbb R^n \f$ becomes + * \f$ \ge \f$. + */ +template +class Multi_critical_filtration { + public: + /** + * @brief Type of the origin of a "lifetime cone". Common with @ref One_critical_filtration "". + */ + using Generator = One_critical_filtration; + using Generators = std::vector; /**< Container type for the filtration values. */ + using iterator = typename Generators::iterator; /**< Iterator type for the generator container. */ + using const_iterator = typename Generators::const_iterator; /**< Const iterator type for the generator container. */ + + // CONSTRUCTORS + + /** + * @brief Default constructor. The constructed value will be either at infinity if `co` is true or at minus infinity + * if `co` is false. + */ + Multi_critical_filtration() : multi_filtration_(_get_default_filtration_value()) {}; + /** + * @brief Constructs a filtration value with one generator and @p n parameters. + * All parameters will be initialized at -inf if `co` is false and at inf if `co` is true. + * + * @warning The generator `{-inf, -inf, ...}`/`{inf, inf, ...}` with \f$ n > 1 \f$ entries is not considered as + * "(minus) infinity" (the resp. methods @ref is_minus_inf() and @ref is_plus_inf() "", as well as the ones of the + * generator, will not return true). The `-inf/inf` are just meant as placeholders, at least one entry should be + * modified by the user. + * Otherwise, either use the static methods @ref minus_inf() or @ref inf(), or set @p n to 1 instead. + * + * @param n Number of parameters. + */ + Multi_critical_filtration(int n) : multi_filtration_(1, Generator(n, _get_default_value())) {}; + /** + * @brief Constructs a filtration value with one generator and @p n parameters. All parameters will be initialized + * with @p value. + * + * @warning If @p value is `inf`, `-inf`, or `NaN`, the generator `{value, value, ...}` with \f$ n > 1 \f$ entries + * is not wrong but will not be considered as respectively "infinity", "minus infinity" or "NaN" (the corresponding + * methods @ref is_plus_inf(), @ref is_minus_inf() and @ref is_nan() will return false). For this purpose, please use + * the static methods @ref inf(), @ref minus_inf() and @ref nan() instead. + * + * @param n Number of parameters. + * @param value Value which will be used for each entry. + */ + Multi_critical_filtration(int n, T value) : multi_filtration_(1, Generator(n, value)) {}; + /** + * @brief Constructs a filtration value with one generator which will be initialzed by the given initializer list. + * + * @param init Initializer list with values for each parameter. + */ + Multi_critical_filtration(std::initializer_list init) : multi_filtration_(1, Generator{init}) {}; + /** + * @brief Constructs a filtration value with one generator which will be initialzed by the given vector. + * + * @param v Vector with values for each parameter. + */ + Multi_critical_filtration(const std::vector &v) : multi_filtration_(1, Generator{v}) {}; + /** + * @brief Constructs a filtration value with one generator to which the given vector is moved to. + * + * @param v Vector with values for each parameter. + */ + Multi_critical_filtration(std::vector &&v) : multi_filtration_(1, Generator{std::move(v)}) {}; + /** + * @brief Constructs filtration value with as many generators than elements in the given vector and initialize + * them with them. + * If the vector is empty, then the filtration value is either initialized at infinity if `co` is true or at + * minus infinity if `co` is false. + * @pre All generators in the vector have to have the same number of parameters, i.e., size. + * Furthermore, the generators have to be a minimal generating set. + * + * @warning If the set of generators is not minimal or not sorted, the behaviour of most methods is undefined. + * It is possible to call @ref simplify() after construction if there is a doubt to ensure this property. + * + * @param v Vector of generators. + */ + Multi_critical_filtration(const std::vector &v) + : multi_filtration_(v.empty() ? _get_default_filtration_value() : v) {}; + /** + * @brief Constructs filtration value with as many generators than elements in the given vector and moves those + * elements to initialize the generators. + * If the vector is empty, then the filtration value is either initialized at infinity if `co` is true or at + * minus infinity if `co` is false. + * @pre All generators in the vector have to have the same number of parameters, i.e., size. + * Furthermore, the generators have to be a minimal generating set. + * + * @warning If the set of generators is not minimal or not sorted, the behaviour of most methods is undefined. + * It is possible to call @ref simplify() after construction if there is a doubt to ensure this property. + * + * @param v Vector of generators. + */ + Multi_critical_filtration(std::vector &&v) + : multi_filtration_(v.empty() ? _get_default_filtration_value() : std::move(v)) {}; + /** + * @brief Constructs a filtration value with one generator initialzed by the range given by the begin and end + * iterators. + * + * @param it_begin Start of the range. + * @param it_end End of the range. + */ + Multi_critical_filtration(typename std::vector::iterator it_begin, typename std::vector::iterator it_end) + : multi_filtration_(Generators(1, {it_begin, it_end})) {}; + /** + * @brief Constructs a filtration value with one generator initialzed by the range given by the begin and end + * const iterators. + * + * @param it_begin Start of the range. + * @param it_end End of the range. + */ + Multi_critical_filtration(typename std::vector::const_iterator it_begin, + typename std::vector::const_iterator it_end) + : multi_filtration_(Generator(1, {it_begin, it_end})) {}; + + // VECTOR-LIKE + + using value_type = T; /**< Entry type. */ + + /** + * @brief Standard operator[]. + */ + Generator &operator[](std::size_t i) { return multi_filtration_[i]; } + /** + * @brief Standard operator[] const. + */ + const Generator &operator[](std::size_t i) const { return multi_filtration_[i]; } + + /** + * @brief Returns begin iterator of the generator range. + * + * @warning If the generator is modified and the new set of generators is not minimal or not sorted, the behaviour + * of most methods is undefined. It is possible to call @ref simplify() after construction if there is a doubt to + * ensure this property. + */ + iterator begin() { return multi_filtration_.begin(); } + + /** + * @brief Returns end iterator of the generator range. + * + * @warning If the generator is modified and the new set of generators is not minimal or not sorted, the behaviour + * of most methods is undefined. It is possible to call @ref simplify() after construction if there is a doubt to + * ensure this property. + */ + iterator end() { return multi_filtration_.end(); } + + /** + * @brief Returns begin const iterator of the generator range. + */ + const_iterator begin() const { return multi_filtration_.begin(); } + + /** + * @brief Returns end const iterator of the generator range. + */ + const_iterator end() const { return multi_filtration_.end(); } + + /** + * @brief Reserves space for the given number of generators in the underlying container. + * + * @param n Number of generators. + */ + void reserve(std::size_t n) { multi_filtration_.reserve(n); } + + // CONVERTERS + + /** + * @brief Casts the object into the type of a generator. + * @pre The filtration value is 1-critical. If there are more than one generator, only the first will be preserved + * and if there is no generator, the method will segfault. + */ + operator Generator() { + if (num_generators() != 1) + throw std::logic_error("Casting a " + std::to_string(num_generators()) + + "-critical filtration value into an 1-critical filtration value."); + return multi_filtration_[0]; + } + + /** + * @brief Casts the object into the type of a generator. + * @pre The filtration value is 1-critical. If there are more than one generator, only the first will be preserved + * and if there is no generator, the method will segfault. + */ + operator const Generator() const { + GUDHI_CHECK(num_generators() == 1, "Casting a " + std::to_string(num_generators()) + + "-critical filtration value into an 1-critical filtration value."); + return multi_filtration_[0]; + } + + // like numpy + /** + * @brief Returns a copy with entries casted into the type given as template parameter. + * + * @tparam U New type for the entries. + * @return Copy with new entry type. + */ + template + Multi_critical_filtration as_type() const { + std::vector> out(num_generators()); + for (std::size_t i = 0u; i < num_generators(); ++i) { + out[i] = multi_filtration_[i].template as_type(); + } + return Multi_critical_filtration(std::move(out)); + } + + // ACCESS + + /** + * @brief Returns a reference to the underlying container storing the generators. + * + * @warning If a generator is modified and the new set of generators is not minimal or not sorted, the behaviour + * of most methods is undefined. It is possible to call @ref simplify() after construction if there is a doubt to + * ensure this property. + */ + const Generators &get_underlying_container() const { return multi_filtration_; } + + /** + * @brief Returns the number of parameters. + */ + std::size_t num_parameters() const { return multi_filtration_[0].num_parameters(); } + + /** + * @brief Returns the number of generators. + */ + std::size_t num_generators() const { return multi_filtration_.size(); } + + /** + * @brief Returns a filtration value for which @ref is_plus_inf() returns `true`. + * + * @return Infinity. + */ + constexpr static Multi_critical_filtration inf() { return Multi_critical_filtration(Generator::inf()); } + + /** + * @brief Returns a filtration value for which @ref is_minus_inf() returns `true`. + * + * @return Minus infinity. + */ + constexpr static Multi_critical_filtration minus_inf() { return Multi_critical_filtration(Generator::minus_inf()); } + + /** + * @brief Returns a filtration value for which @ref is_nan() returns `true`. + * + * @return NaN. + */ + constexpr static Multi_critical_filtration nan() { return Multi_critical_filtration(Generator::nan()); } + + // DESCRIPTORS + + // TODO: Accept {{-inf, -inf, ...},...} / {{inf, inf, ...},...} / {{NaN, NaN, ...},...} as resp. -inf / inf / NaN. + + /** + * @brief Returns `true` if and only if the filtration value is considered as infinity. + */ + bool is_plus_inf() const { return multi_filtration_.size() == 1 && multi_filtration_[0].is_plus_inf(); } + + /** + * @brief Returns `true` if and only if the filtration value is considered as minus infinity. + */ + bool is_minus_inf() const { return multi_filtration_.size() == 1 && multi_filtration_[0].is_minus_inf(); } + + /** + * @brief Returns `true` if and only if the filtration value is considered as NaN. + */ + bool is_nan() const { return multi_filtration_.size() == 1 && multi_filtration_[0].is_nan(); } + + /** + * @brief Returns `true` if and only if the filtration value is non-empty and is not considered as infinity, + * minus infinity or NaN. + */ + bool is_finite() const { + if (multi_filtration_.size() > 1) return true; + return multi_filtration_[0].is_finite(); + } + + // COMPARAISON OPERATORS + + // TODO : this costs a lot... optimize / cheat in some way for python ? + + /** + * @brief Returns `true` if and only if the positive cones generated by @p b are strictly contained in the + * positive cones generated by @p a. + * If @p a and @p b are both not infinite or NaN, they have to have the same number of parameters. + * + * Note that not all filtration values are comparable. That is, \f$ a > b \f$ and \f$ b > a \f$ returning both false + * does **not** imply \f$ a == b \f$. + */ + friend bool operator<(const Multi_critical_filtration &a, const Multi_critical_filtration &b) { + for (std::size_t i = 0u; i < b.multi_filtration_.size(); ++i) { + // for each generator in b, verify if it is strictly in the cone of at least one generator of a + bool isContained = false; + for (std::size_t j = 0u; j < a.multi_filtration_.size() && !isContained; ++j) { + // lexicographical order, so if a[j][0] dom b[j][0], than a[j'] can never strictly contain b[i] for all j' > j. + if (_first_dominates(a.multi_filtration_[j], b.multi_filtration_[i])) return false; + isContained = _strictly_contains(a.multi_filtration_[j], b.multi_filtration_[i]); + } + if (!isContained) return false; + } + return true; + } + + /** + * @brief Returns `true` if and only if the positive cones generated by @p a are strictly contained in the + * positive cones generated by @p b. + * If @p a and @p b are both not infinite or NaN, they have to have the same number of parameters. + * + * Note that not all filtration values are comparable. That is, \f$ a > b \f$ and \f$ b > a \f$ returning both false + * does **not** imply \f$ a == b \f$. + */ + friend bool operator>(const Multi_critical_filtration &a, const Multi_critical_filtration &b) { return b < a; } + + /** + * @brief Returns `true` if and only if the positive cones generated by @p b are contained in or are (partially) + * equal to the positive cones generated by @p a. + * If @p a and @p b are both not infinite or NaN, they have to have the same number of parameters. + * + * Note that not all filtration values are comparable. That is, \f$ a \le b \f$ and \f$ b \le a \f$ can both return + * `false`. + */ + friend bool operator<=(const Multi_critical_filtration &a, const Multi_critical_filtration &b) { + // check if this curves is below other's curve + // ie for each guy in this, check if there is a guy in other that dominates him + for (std::size_t i = 0u; i < b.multi_filtration_.size(); ++i) { + // for each generator in b, verify if it is in the cone of at least one generator of a + bool isContained = false; + for (std::size_t j = 0u; j < a.multi_filtration_.size() && !isContained; ++j) { + // lexicographical order, so if a[j][0] strictly dom b[j][0], than a[j'] can never contain b[i] for all j' > j. + if (_first_strictly_dominates(a.multi_filtration_[j], b.multi_filtration_[i])) return false; + isContained = _contains(a.multi_filtration_[j], b.multi_filtration_[i]); + } + if (!isContained) return false; + } + return true; + } + + /** + * @brief Returns `true` if and only if the positive cones generated by @p a are contained in or are (partially) + * equal to the positive cones generated by @p b. + * If @p a and @p b are both not infinite or NaN, they have to have the same number of parameters. + * + * Note that not all filtration values are comparable. That is, \f$ a \ge b \f$ and \f$ b \ge a \f$ can both return + * `false`. + */ + friend bool operator>=(const Multi_critical_filtration &a, const Multi_critical_filtration &b) { return b <= a; } + + /** + * @brief Returns `true` if and only if for each \f$ i \f$, \f$ a[i] \f$ is equal to \f$ b[i] \f$. + */ + friend bool operator==(const Multi_critical_filtration &a, const Multi_critical_filtration &b) { + // assumes lexicographical order for both + return a.multi_filtration_ == b.multi_filtration_; + } + + /** + * @brief Returns `true` if and only if \f$ a == b \f$ returns `false`. + */ + friend bool operator!=(const Multi_critical_filtration &a, const Multi_critical_filtration &b) { return !(a == b); } + + // MODIFIERS + + /** + * @brief Sets the number of generators. If there were less generators before, new empty generators are constructed. + * If there were more generators before, the exceed of generators is destroyed (any generator with index higher or + * equal than @p n to be more precise). If @p n is zero, the methods does nothing. A filtration value should never + * be empty. + * + * @warning All empty generators have 0 parameters. This can be problematic for some methods if there are also + * non empty generators in the container. Make sure to fill them with real generators or to remove them before + * using those methods. + * + * @warning Be sure to call @ref simplify if necessary after setting all the generators. Most methods will have an + * undefined behaviour if the set of generators is not minimal or sorted. + * + * @param n New number of generators. + */ + void set_num_generators(std::size_t n) { + if (n == 0) return; + multi_filtration_.resize(n); + } + + /** + * @brief Sets all generators to the least common upper bound between the current generator value and the given value. + * + * More formally, it pushes the current filtration value to the cone \f$ \{ y \in \mathbb R^n : y \ge x \} \f$ + * originating in \f$ x \f$. The resulting values corresponds to the generators of the intersection of this cone + * with the union of positive cones generated by the old generators. + * + * @param x The target filtration value towards which to push. + */ + void push_to_least_common_upper_bound(const Generator &x) { + if (this->is_plus_inf() || this->is_nan() || x.is_nan() || x.is_minus_inf()) return; + + GUDHI_CHECK(x.is_plus_inf() || x.num_parameters() == multi_filtration_[0].num_parameters() || !is_finite(), + "Pushing to a filtration value with different number of parameters."); + + if (x.is_plus_inf() || this->is_minus_inf()) { + multi_filtration_ = {x}; + return; + } + for (auto &g : *this) { + g.push_to_least_common_upper_bound(x); + } + + simplify(); + } + + /** + * @brief Sets all generators to the greatest common lower bound between the current generator value and the given + * value. + * + * More formally, it pulls the current filtration value to the cone \f$ \{ y \in \mathbb R^n : y \le x \} \f$ + * originating in \f$ x \f$. The resulting values corresponds to the generators of the intersection of this cone + * with the union of negative cones generated by the old generators. + * + * @param x The target filtration value towards which to pull. + */ + void pull_to_greatest_common_lower_bound(const Generator &x) { + if (x.is_plus_inf() || this->is_nan() || x.is_nan() || this->is_minus_inf()) return; + + GUDHI_CHECK(x.is_minus_inf() || x.num_parameters() == multi_filtration_[0].num_parameters() || !is_finite(), + "Pulling to a filtration value with different number of parameters."); + + if (this->is_plus_inf() || x.is_minus_inf()) { + multi_filtration_ = {x}; + return; + } + for (auto &g : *this) { + g.pull_to_greatest_common_lower_bound(x); + } + + simplify(); + } + + /** + * @brief Adds the given generator to the filtration value such that the sets remains minimal. + * It is therefore possible that the generator is ignored if it does not generated any new lifetime or that + * old generators disappear if they are overshadowed by the new one. + * @pre If all are finite, the new generator has to have the same number of parameters than the others. + * + * @param x New generator to add. + * @return true If and only if the generator is actually added to the set of generators. + * @return false Otherwise. + */ + bool add_generator(const Generator &x) { + GUDHI_CHECK(x.num_parameters() == multi_filtration_[0].num_parameters() || !is_finite() || !x.is_finite(), + "Cannot add a generator with different number of parameters."); + + std::size_t end = multi_filtration_.size(); + + if (_generator_can_be_added(x, 0, end)) { + multi_filtration_.resize(end); + multi_filtration_.push_back(x); + std::sort(multi_filtration_.begin(), multi_filtration_.end(), Is_strictly_smaller_lexicographically()); + return true; + } + + return false; + } + + /** + * @brief Adds the given generator to the filtration value without any verifications or simplifications. + * + * @warning If the resulting set of generators is not minimal after modification, some methods will have an + * undefined behaviour. Be sure to call @ref simplify() before using them. + * + * @param x + */ + void add_guaranteed_generator(const Generator &x) { multi_filtration_.push_back(x); } + + /* + * Same as `compute_coordinates_in_grid`, but does the operation in-place. + */ + + /** + * @brief Projects the filtration value into the given grid. If @p coordinate is false, the entries are set to + * the nearest upper bound value with the same parameter in the grid and the new generators are simplified and + * ordered. Otherwise, the entries are set to the indices of those nearest upper bound values. In this case, + * no simplification or sort are done, such that the new coordinates have a one by one correspondence with the + * positions of the old generators. + * The grid has to be represented as a vector of ordered ranges of values convertible into `T`. An index + * \f$ i \f$ of the vector corresponds to the same parameter as the index \f$ i \f$ in a generator. + * The ranges correspond to the possible values of the parameters, ordered by increasing value, forming therefore + * all together a 2D grid. + * + * @tparam oned_array A range of values convertible into `T` ordered by increasing value. Has to implement + * a begin, end and operator[] method. + * @param grid Vector of @p oned_array with size at least number of filtration parameters. + * @param coordinate If true, the values are set to the coordinates of the projection in the grid. If false, + * the values are set to the values at the coordinates of the projection. + */ + template + void project_onto_grid(const std::vector &grid, bool coordinate = true) { + GUDHI_CHECK(grid.size() >= num_parameters(), + "The grid should not be smaller than the number of parameters in the filtration value."); + + for (auto &x : multi_filtration_) { + x.project_onto_grid(grid, coordinate); + } + + if (!coordinate) simplify(); + } + + /** + * @brief Removes all empty generators from the filtration value. If @p include_infinities is true, it also + * removes the generators at infinity or minus infinity. + * If the set of generators is empty after removals, it is set to minus infinity if `co` is false or to infinity + * if `co` is true. + * + * @param include_infinities If true, removes also infinity values. + */ + void remove_empty_generators(bool include_infinities = false) { + multi_filtration_.erase(std::remove_if(multi_filtration_.begin(), multi_filtration_.end(), + [include_infinities](const Generator &a) { + return a.empty() || + ((include_infinities) && (a.is_plus_inf() || a.is_minus_inf())); + }), + multi_filtration_.end()); + std::sort(multi_filtration_.begin(), multi_filtration_.end(), Is_strictly_smaller_lexicographically()); + if (multi_filtration_.empty()) multi_filtration_.push_back(Generator{_get_default_value()}); + } + + /** + * @brief Simplifies the current set of generators such that it becomes minimal. Also orders it in increasing + * lexicographical order. Only necessary if generators were added "by hand" without verification either trough the + * constructor or with @ref add_guaranteed_generator "", etc. + */ + void simplify() { + std::size_t end = 0; + + for (std::size_t curr = 0; curr < multi_filtration_.size(); ++curr) { + if (_generator_can_be_added(multi_filtration_[curr], 0, end)) { + std::swap(multi_filtration_[end], multi_filtration_[curr]); + ++end; + } + } + + multi_filtration_.resize(end); + std::sort(multi_filtration_.begin(), multi_filtration_.end(), Is_strictly_smaller_lexicographically()); + } + + // FONCTIONNALITIES + + /** + * @brief Returns a generator with the minimal values of all parameters in any generator of the given filtration + * value. That is, the greatest lower bound of all generators. + */ + friend Generator factorize_below(const Multi_critical_filtration &f) { + if (f.num_generators() == 0) return Generator(); + Generator result(f.num_parameters(), Generator::T_inf); + for (const auto &g : f) { + if (g.is_nan() || g.is_minus_inf()) return g; + if (g.is_plus_inf()) continue; + for (std::size_t i = 0; i < f.num_parameters(); ++i) { + result[i] = std::min(result[i], g[i]); + } + } + return result; + } + + /** + * @brief Returns a generator with the maximal values of all parameters in any generator of the given filtration + * value. That is, the least upper bound of all generators. + */ + friend Generator factorize_above(const Multi_critical_filtration &f) { + if (f.num_generators() == 0) return Generator(); + Generator result(f.num_parameters(), -Generator::T_inf); + for (auto &g : f) { + if (g.is_nan() || g.is_plus_inf()) return g; + if (g.is_minus_inf()) continue; + for (std::size_t i = 0; i < g.num_parameters(); ++i) { + result[i] = std::max(result[i], g[i]); + } + } + return result; + } + + /** + * @brief Computes the smallest (resp. the greatest if `co` is true) scalar product of the all generators with the + * given vector. + * + * @tparam U Arithmetic type of the result. Default value: `T`. + * @param f Filtration value. + * @param x Vector of coefficients. + * @return Scalar product of @p f with @p x. + */ + template + friend U compute_linear_projection(const Multi_critical_filtration &f, const std::vector &x) { + if constexpr (co) { + U projection = std::numeric_limits::lowest(); + for (const auto &y : f) { + projection = std::max(projection, compute_linear_projection(y, x)); + } + return projection; + } else { + U projection = std::numeric_limits::max(); + for (const auto &y : f) { + projection = std::min(projection, compute_linear_projection(y, x)); + } + return projection; + } + } + + /** + * @brief Computes the coordinates in the given grid, corresponding to the nearest upper bounds of the entries + * in the given filtration value. + * The grid has to be represented as a vector of vectors of ordered values convertible into `out_type`. An index + * \f$ i \f$ of the vector corresponds to the same parameter as the index \f$ i \f$ in a generator. + * The inner vectors correspond to the possible values of the parameters, ordered by increasing value, + * forming therefore all together a 2D grid. + * + * @tparam out_type Signed arithmetic type. Default value: std::int32_t. + * @tparam U Type which is convertible into `out_type`. + * @param f Filtration value to project. + * @param grid Vector of vectors to project into. + * @return Filtration value \f$ out \f$ whose entry correspond to the indices of the projected values. That is, + * the projection of \f$ f[i] \f$ is \f$ grid[i][out[i]] \f$ before simplification (if two generators were + * projected to the same point, the doubles are removed in the output). + */ + template + friend Multi_critical_filtration compute_coordinates_in_grid(Multi_critical_filtration f, + const std::vector> &grid) { + // TODO: by replicating the code of the 1-critical "project_onto_grid", this could be done with just one copy + // instead of two. But it is not clear if it is really worth it, i.e., how much the change in type is really + // necessary in the use cases. To see later. + f.project_onto_grid(grid); + if constexpr (std::is_same_v) { + return f; + } else { + return f.as_type(); + } + } + + /** + * @brief Computes the values in the given grid corresponding to the coordinates given by the given filtration + * value. That is, if \f$ out \f$ is the result, \f$ out[i] = grid[i][f[i]] \f$. Assumes therefore, that the + * values stored in the filtration value corresponds to indices existing in the given grid. + * + * @tparam U Signed arithmetic type. + * @param f Filtration value storing coordinates compatible with `grid`. + * @param grid Vector of vector. + * @return Filtration value \f$ out \f$ whose entry correspond to \f$ out[i] = grid[i][f[i]] \f$ before + * simplification (the output is simplified). + */ + template + friend Multi_critical_filtration evaluate_coordinates_in_grid(const Multi_critical_filtration &f, + const std::vector> &grid) { + Multi_critical_filtration out; + out.set_num_generators(f.num_generators()); + for (std::size_t i = 0; i < f.num_generators(); ++i) { + out[i] = evaluate_coordinates_in_grid(f[i], grid); + } + out.simplify(); + return out; + } + + // UTILITIES + + /** + * @brief Outstream operator. + */ + friend std::ostream &operator<<(std::ostream &stream, const Multi_critical_filtration &f) { + if (f.is_plus_inf()) { + stream << "[inf, ..., inf]"; + return stream; + } + if (f.is_minus_inf()) { + stream << "[-inf, ..., -inf]"; + return stream; + } + if (f.is_nan()) { + stream << "[NaN]"; + return stream; + } + stream << "(k=" << f.multi_filtration_.size() << ")["; + for (const auto &val : f) { + stream << val << "; "; + } + if (f.multi_filtration_.size() > 0) { + stream << "\b" + << "\b"; + } + stream << "]"; + return stream; + } + + public: + /** + * @brief Indicates if the class manages multi-critical filtration values. + */ + constexpr static const bool is_multi_critical = true; + + private: + Generators multi_filtration_; /**< Container for generators. */ + + struct Is_strictly_smaller_lexicographically { + //assumes both generators have the same length if not infinite/nan. + bool operator()(const Generator &g1, const Generator &g2) { + if (g1.is_nan() || g2.is_nan()) return !g1.is_nan(); + if (g1.is_plus_inf()) return false; + if (g2.is_plus_inf()) return true; + if (g1.is_minus_inf()) return false; + if (g2.is_minus_inf()) return true; + + // g1 and g2 have to finite and of the same size + for (std::size_t i = 0; i < g1.size(); ++i) { + if (g1[i] != g2[i]) return g1[i] < g2[i]; + } + return false; + } + }; + + constexpr static T _get_default_value() { return co ? Generator::T_inf : -Generator::T_inf; } + + constexpr static Generators _get_default_filtration_value() { return Generators{Generator{_get_default_value()}}; } + + /** + * @brief Verifies if @p b is strictly contained in the positive cone originating in `a`. + */ + static bool _strictly_contains(const Generator &a, const Generator &b) { + if constexpr (co) + return a > b; + else { + return a < b; + } + } + /** + * @brief Verifies if @p b is contained in the positive cone originating in `a`. + */ + static bool _contains(const Generator &a, const Generator &b) { + if constexpr (co) + return a >= b; + else { + return a <= b; + } + } + + static bool _first_strictly_dominates(const Generator& a, const Generator& b){ + if constexpr (co){ + return !a.empty() && !b.empty() && a[0] < b[0]; + } else { + return !a.empty() && !b.empty() && a[0] > b[0]; + } + } + + static bool _first_dominates(const Generator& a, const Generator& b){ + if constexpr (co){ + return !a.empty() && !b.empty() && a[0] <= b[0]; + } else { + return !a.empty() && !b.empty() && a[0] >= b[0]; + } + } + + enum class Rel { EQUAL, DOMINATES, IS_DOMINATED, NONE }; + + static Rel _get_domination_relation(const Generator &a, const Generator &b) { + if (a.is_nan() || b.is_nan()) return Rel::NONE; + + GUDHI_CHECK(a.size() == b.size(), + "Two generators in the same k-critical value have to have the same numbers of parameters."); + + bool equal = true; + bool allGreater = true; + bool allSmaller = true; + for (unsigned int i = 0; i < a.size(); ++i) { + if (a[i] < b[i]) { + if (!allSmaller) return Rel::NONE; + equal = false; + allGreater = false; + } else if (a[i] > b[i]) { + if (!allGreater) return Rel::NONE; + equal = false; + allSmaller = false; + } + } + if (equal) return Rel::EQUAL; + + if constexpr (co) { + if (allSmaller) return Rel::DOMINATES; + return Rel::IS_DOMINATED; + } else { + if (allGreater) return Rel::DOMINATES; + return Rel::IS_DOMINATED; + } + } + + // assumes between 'curr' and 'end' everything is simplified: + // no nan values and if there is an inf/-inf, then 'end - curr == 1' + // modifies multi_filtration_ only if true is returned. + bool _generator_can_be_added(const Generator &x, std::size_t curr, std::size_t &end) { + if (x.empty() || x.is_nan()) return false; + + // assumes that everything between curr and end is simplified + // so, only multi_filtration_[curr] can be at inf or -inf. + if constexpr (co) { + if (multi_filtration_[curr].is_plus_inf() || (x.is_minus_inf() && end - curr != 0)) { + return false; + } + if (multi_filtration_[curr].is_minus_inf()) { + if (x.is_minus_inf()) { + return false; + } + end = curr; + return true; + } + if (x.is_plus_inf()) { + if (multi_filtration_[curr].is_plus_inf()) return false; + end = curr; + return true; + } + } else { + if (multi_filtration_[curr].is_minus_inf() || (x.is_plus_inf() && end - curr != 0)) { + return false; + } + if (multi_filtration_[curr].is_plus_inf()) { + if (x.is_plus_inf()) { + return false; + } + end = curr; + return true; + } + if (x.is_minus_inf()) { + if (multi_filtration_[curr].is_minus_inf()) return false; + end = curr; + return true; + } + } + + while (curr != end) { + Rel res = _get_domination_relation(multi_filtration_[curr], x); + if (res == Rel::IS_DOMINATED || res == Rel::EQUAL) return false; // x dominates or is equal + if (res == Rel::DOMINATES) { // x is dominated + --end; + std::swap(multi_filtration_[curr], multi_filtration_[end]); + } else { // no relation + ++curr; + } + } + return true; + } +}; + +} // namespace Gudhi::multi_filtration + +namespace std { + +template +class numeric_limits> { + public: + static constexpr bool has_infinity = true; + + static constexpr Gudhi::multi_filtration::Multi_critical_filtration infinity() noexcept { + return Gudhi::multi_filtration::Multi_critical_filtration::inf(); + }; + + // non-standard + static constexpr Gudhi::multi_filtration::Multi_critical_filtration minus_infinity() noexcept { + return Gudhi::multi_filtration::Multi_critical_filtration::minus_inf(); + }; + + static constexpr Gudhi::multi_filtration::Multi_critical_filtration max() noexcept(false) { + throw std::logic_error( + "The maximal value cannot be represented with no finite numbers of generators." + "Use `max(number_of_generators, number_of_parameters)` instead"); + }; + + // non-standard, so I don't want to define default values. + static constexpr Gudhi::multi_filtration::Multi_critical_filtration max(unsigned int g, unsigned int n) noexcept { + std::vector::Generator> v( + g, std::vector(n, std::numeric_limits::max())); + return Gudhi::multi_filtration::Multi_critical_filtration(std::move(v)); + }; + + static constexpr Gudhi::multi_filtration::Multi_critical_filtration quiet_NaN() noexcept { + return Gudhi::multi_filtration::Multi_critical_filtration::nan(); + }; +}; + +} // namespace std + +#endif // MULTI_CRITICAL_FILTRATIONS_H_ diff --git a/src/Multi_filtration/include/gudhi/One_critical_filtration.h b/src/Multi_filtration/include/gudhi/One_critical_filtration.h new file mode 100644 index 0000000000..8b44fbfcb0 --- /dev/null +++ b/src/Multi_filtration/include/gudhi/One_critical_filtration.h @@ -0,0 +1,1378 @@ +/* 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): David Loiseaux + * + * Copyright (C) 2023 Inria + * + * Modification(s): + * - 2024/08 Hannah Schreiber: Generalization to all signed arithmetic types for T + doc + * - YYYY/MM Author: Description of the modification + */ + +/** + * @file One_critical_filtration.h + * @author David Loiseaux + * @brief Contains the @ref Gudhi::multi_filtration::One_critical_filtration class. + */ + +#ifndef ONE_CRITICAL_FILTRATIONS_H_ +#define ONE_CRITICAL_FILTRATIONS_H_ + +#include //std::lower_bound +#include //std::isnan +#include //std::size_t +#include //std::int32_t +#include //std::ostream +#include //std::numerical_limits +#include //std::logic_error +#include + +#include + +namespace Gudhi::multi_filtration { + +/** + * @class One_critical_filtration one_critical_filtration.h gudhi/one_critical_filtration.h + * @ingroup multi_filtration + * + * @brief Class encoding the apparition time, i.e., filtration value of an object + * (e.g., simplex, cell, abstract algebraic generator) in the setting of 1-critical multiparameter filtrations. + * The class can be used as a vector whose indices correspond to one parameter each. + * It also follows numpy-like broadcast semantic. + * + * @details Inherits of `std::vector`. Overloads `std::numeric_limits` such that: + * - `std::numeric_limits >::has_infinity` returns `true`, + * - `std::numeric_limits >::infinity()` returns @ref One_critical_filtration::inf() "", + * - `std::numeric_limits >::minus_infinity()` returns + * @ref One_critical_filtration::minus_inf() "", + * - `std::numeric_limits >::max()` throws, + * - `std::numeric_limits >::max(n)` returns a @ref One_critical_filtration with `n` + * parameters evaluated at value `std::numeric_limits::max()`, + * - `std::numeric_limits >::quiet_NaN()` returns @ref One_critical_filtration::nan() "". + * + * One critical simplicial filtrations are filtrations such that the lifetime of each object is a positive cone, e.g. + * - \f$ \{ x \in \mathbb R^2 : x>=(1,2)\} \f$ is valid, while + * - \f$ \{ x \in \mathbb R^2 : x>=(1,2)\} \cap \{x \in \mathbb R^2 : x>=(2,1)\} \f$ is not. + * + * If the lifetime corresponds to a union of such positive cones, the filtration is called a multi-critical filtration. + * For those cases, use @ref Multi_critical_filtration instead. + * + * @tparam T Arithmetic type of an entry for one parameter of the filtration value. Has to be **signed** and + * to implement `std::isnan(T)`, `std::numeric_limits::has_quiet_NaN`, `std::numeric_limits::quiet_NaN()`, + * `std::numeric_limits::has_infinity`, `std::numeric_limits::infinity()` and `std::numeric_limits::max()`. + * If `std::numeric_limits::has_infinity` returns `false`, a call to `std::numeric_limits::infinity()` + * can simply throw. Examples are the native types `double`, `float` and `int`. + */ +template +class One_critical_filtration : public std::vector { + private: + using Base = std::vector; + + public: + /** + * @brief Type of the origin of a "lifetime cone", i.e., of a one-critical filtration value. + * Common with @ref Multi_critical_filtration "". In the 1-critical case, simply the class it-self. + * + * @tparam U Type of a value for one parameter within the filtration value. Default value: `T`. + */ + template + using Generator = One_critical_filtration; + + // CONSTRUCTORS + + /** + * @brief Default constructor. Constructs a value at minus infinity. + */ + One_critical_filtration() : Base{-T_inf} {}; + /** + * @brief Constructs a vector of the size of the given number of parameters with -inf as value for each entry. + * + * @warning The vector `{-inf, -inf, ...}` with \f$ n > 1 \f$ entries is not considered as "minus infinity" (the + * method @ref is_minus_inf() will not return true). The `-inf` are just meant as placeholders, at least one entry + * should be modified by the user. Otherwise, either use the static method @ref minus_inf() or set @p n to 1 instead. + * + * @param n Number of parameters. + */ + One_critical_filtration(int n) : Base(n, -T_inf) {}; + /** + * @brief Constructs a vector of the size of the given number of parameters and the given value for each entry. + * + * @warning If @p value is `inf`, `-inf`, or `NaN`, the vector `{value, value, ...}` with \f$ n > 1 \f$ entries + * is not wrong but will not be considered as respectively "infinity", "minus infinity" or "NaN" (the corresponding + * methods @ref is_plus_inf(), @ref is_minus_inf() and @ref is_nan() will return false). For this purpose, please use + * the static methods @ref inf(), @ref minus_inf() and @ref nan() instead. + * + * @param n Number of parameters. + * @param value Value which will be used for each entry. + */ + One_critical_filtration(int n, T value) : Base(n, value) {}; + /** + * @brief Construct a vector from the given initializer list. + * + * @param init Initializer list with values for each parameter. + */ + One_critical_filtration(std::initializer_list init) : Base(init) {}; + /** + * @brief Construct a vector from the given vector. + * + * @param v Vector with values for each parameter. + */ + One_critical_filtration(const std::vector &v) : Base(v) {}; + /** + * @brief Construct a vector from the given vector by moving it to the new vector. + * + * @param v Vector with values for each parameter. + */ + One_critical_filtration(std::vector &&v) : Base(std::move(v)) {}; + /** + * @brief Construct a vector from the range given by the begin and end iterators. + * + * @param it_begin Start of the range. + * @param it_end End of the range. + */ + One_critical_filtration(typename std::vector::iterator it_begin, typename std::vector::iterator it_end) + : Base(it_begin, it_end) {}; + /** + * @brief Construct a vector from the range given by the begin and end const iterators. + * + * @param it_begin Start of the range. + * @param it_end End of the range. + */ + One_critical_filtration(typename std::vector::const_iterator it_begin, + typename std::vector::const_iterator it_end) + : Base(it_begin, it_end) {}; + + // HERITAGE + + using std::vector::operator[]; /**< Inheritance of entry access. */ + using value_type = T; /**< Entry type. */ + + // CONVERTERS + + // like numpy + /** + * @brief Returns a copy with entries casted into the type given as template parameter. + * + * @tparam U New type for the entries. + * @return Copy with new entry type. + */ + template + One_critical_filtration as_type() const { + One_critical_filtration out(0); + out.reserve(Base::size()); + for (std::size_t i = 0u; i < Base::size(); i++) out.push_back(static_cast(Base::operator[](i))); + return out; + } + + // ACCESS + + /** + * @brief Returns the number of parameters of the finite filtration value. If the value is "inf", "-inf" or "NaN", + * returns 1. + * + * @return Number of parameters. + */ + std::size_t num_parameters() const { return Base::size(); } + + /** + * @brief Returns a filtration value for which @ref is_plus_inf() returns `true`. + * + * @return Infinity. + */ + constexpr static One_critical_filtration inf() { return {T_inf}; } + + /** + * @brief Returns a filtration value for which @ref is_minus_inf() returns `true`. + * + * @return Minus infinity. + */ + constexpr static One_critical_filtration minus_inf() { return {-T_inf}; } + + /** + * @brief Returns a filtration value for which @ref is_nan() returns `true`. + * + * @return NaN. + */ + constexpr static One_critical_filtration nan() { + if constexpr (std::numeric_limits::has_quiet_NaN) { + return {std::numeric_limits::quiet_NaN()}; + } else { + return One_critical_filtration(0); // to differentiate it from 0, an empty filtration value can't do much anyway. + } + } + + // DESCRIPTORS + + // TODO: Accept {-inf, -inf, ...} / {inf, inf, ...} / {NaN, NaN, ...} as resp. -inf / inf / NaN. + + /** + * @brief Returns `true` if and only if the filtration value is considered as infinity. + */ + bool is_plus_inf() const { + if (Base::size() != 1) return false; + return (Base::operator[](0) == T_inf); + } + + /** + * @brief Returns `true` if and only if the filtration value is considered as minus infinity. + */ + bool is_minus_inf() const { + if (Base::size() != 1) return false; + return (Base::operator[](0) == -T_inf); + } + + /** + * @brief Returns `true` if and only if the filtration value is considered as NaN. + */ + bool is_nan() const { + if constexpr (std::numeric_limits::has_quiet_NaN) { + if (Base::size() != 1) return false; + return is_nan_(Base::operator[](0)); + } else { + return Base::empty(); + } + } + + /** + * @brief Returns `true` if and only if the filtration value is non-empty and is not considered as infinity, + * minus infinity or NaN. + */ + bool is_finite() const { + if (Base::size() > 1) return true; + if (Base::size() == 0) return false; + auto first_value = Base::operator[](0); // TODO : Maybe check all entries ? + if (is_nan_(first_value) || first_value == -T_inf || first_value == T_inf) return false; + return true; + } + + // COMPARAISON OPERATORS + + /** + * @brief Returns `true` if and only if for each \f$ i \f$, \f$ a[i] \f$ is strictly smaller than \f$ b[i] \f$. + * If @p a and @p b are both not infinite or NaN, they have to have the same number of parameters. + * + * Note that not all filtration values are comparable. That is, \f$ a < b \f$ and \f$ b < a \f$ returning both false + * does **not** imply \f$ a == b \f$. + */ + friend bool operator<(const One_critical_filtration &a, const One_critical_filtration &b) { + if (a.is_plus_inf() || a.is_nan() || b.is_nan() || b.is_minus_inf()) return false; + if (b.is_plus_inf() || a.is_minus_inf()) return true; + bool isSame = true; + auto n = a.size(); + GUDHI_CHECK(a.size() == b.size(), "Two filtration values with different number of parameters are not comparable."); + for (auto i = 0u; i < n; ++i) { + if (a[i] > b[i]) return false; + if (isSame && a[i] != b[i]) isSame = false; + } + return !isSame; + } + + /** + * @brief Returns `true` if and only if for each \f$ i \f$, \f$ a[i] \f$ is smaller or equal than \f$ b[i] \f$. + * If @p a and @p b are both not infinite or NaN, they have to have the same number of parameters. + * + * Note that not all filtration values are comparable. That is, \f$ a \le b \f$ and \f$ b \le a \f$ can both return + * `false`. + */ + friend bool operator<=(const One_critical_filtration &a, const One_critical_filtration &b) { + if (a.is_nan() || b.is_nan()) return false; + if (b.is_plus_inf() || a.is_minus_inf()) return true; + if (a.is_plus_inf() || b.is_minus_inf()) return false; + auto n = a.size(); + GUDHI_CHECK(a.size() == b.size(), "Two filtration values with different number of parameters are not comparable."); + for (std::size_t i = 0u; i < n; ++i) { + if (a[i] > b[i]) return false; + } + return true; + } + + /** + * @brief Returns `true` if and only if for each \f$ i \f$, \f$ a[i] \f$ is strictly greater than \f$ b[i] \f$. + * If @p a and @p b are both not infinite or NaN, they have to have the same number of parameters. + * + * Note that not all filtration values are comparable. That is, \f$ a > b \f$ and \f$ b > a \f$ returning both false + * does **not** imply \f$ a == b \f$. + */ + friend bool operator>(const One_critical_filtration &a, const One_critical_filtration &b) { return b < a; } + + /** + * @brief Returns `true` if and only if for each \f$ i \f$, \f$ a[i] \f$ is greater or equal than \f$ b[i] \f$. + * If @p a and @p b are both not infinite or NaN, they have to have the same number of parameters. + * + * Note that not all filtration values are comparable. That is, \f$ a \ge b \f$ and \f$ b \ge a \f$ can both return + * `false`. + */ + friend bool operator>=(const One_critical_filtration &a, const One_critical_filtration &b) { return b <= a; } + + /** + * @brief Returns `true` if and only if for each \f$ i \f$, \f$ a[i] \f$ is equal to \f$ b[i] \f$. + */ + friend bool operator==(const One_critical_filtration &a, const One_critical_filtration &b) { + return static_cast(a) == static_cast(b); + } + + /** + * @brief Returns `true` if and only if \f$ a == b \f$ returns `false`. + */ + friend bool operator!=(const One_critical_filtration &a, const One_critical_filtration &b) { return !(a == b); } + + // ARITHMETIC OPERATORS + + // opposite + /** + * @brief Returns a filtration value such that an entry at index \f$ i \f$ is equal to \f$ -f[i] \f$. + * + * Used conventions: + * - \f$ -NaN = NaN \f$. + * + * @param f Value to opposite. + * @return The opposite of @p f. + */ + friend One_critical_filtration operator-(const One_critical_filtration &f) { + One_critical_filtration result(0); + result.reserve(f.size()); + for (auto val : f) { + result.push_back(-val); + } + return result; + } + + // subtraction + /** + * @brief Returns a filtration value such that an entry at index \f$ i \f$ is equal to \f$ a[i] - b[i] \f$. + * If @p a and @p b are both not infinite or NaN, they have to have the same number of parameters. + * + * Used conventions: + * - \f$ inf - inf = NaN \f$, + * - \f$ -inf - (-inf) = NaN \f$, + * - \f$ NaN - b = NaN \f$, + * - \f$ a - NaN = NaN \f$. + * + * If `std::numeric_limits::has_quiet_NaN` is false, then the returned filtration value will be @ref nan() + * if any operation results in NaN, not only if all operations result in NaN. + * + * @param a First element of the subtraction. + * @param b Second element of the subtraction. + * @return Entry-wise \f$ a - b \f$. + */ + friend One_critical_filtration operator-(One_critical_filtration a, const One_critical_filtration &b) { + a -= b; + return a; + } + + /** + * @brief Returns a filtration value such that an entry at index \f$ i \f$ is equal to \f$ f[i] - val \f$. + * + * Used conventions: + * - \f$ inf - inf = NaN \f$, + * - \f$ -inf - (-inf) = NaN \f$, + * - \f$ NaN - b = NaN \f$, + * - \f$ a - NaN = NaN \f$. + * + * If `std::numeric_limits::has_quiet_NaN` is false, then the returned filtration value will be @ref nan() + * if any operation results in NaN, not only if all operations result in NaN. + * + * @param f First element of the subtraction. + * @param val Second element of the subtraction. + * @return Entry-wise \f$ f - val \f$. + */ + friend One_critical_filtration operator-(One_critical_filtration f, const T &val) { + f -= val; + return f; + } + + /** + * @brief Returns a filtration value such that an entry at index \f$ i \f$ is equal to \f$ val - f[i] \f$. + * + * Used conventions: + * - \f$ inf - inf = NaN \f$, + * - \f$ -inf - (-inf) = NaN \f$, + * - \f$ NaN - b = NaN \f$, + * - \f$ a - NaN = NaN \f$. + * + * If `std::numeric_limits::has_quiet_NaN` is false, then the returned filtration value will be @ref nan() + * if any operation results in NaN, not only if all operations result in NaN. + * + * @param val First element of the subtraction. + * @param f Second element of the subtraction. + * @return Entry-wise \f$ val - f \f$. + */ + friend One_critical_filtration operator-(const T &val, One_critical_filtration f) { + // TODO: in one go + f = -f; + f += val; + return f; + } + + /** + * @brief Modifies the first parameters such that an entry at index \f$ i \f$ is equal to + * \f$ result[i] - to\_subtract[i] \f$. + * If @p result and @p to_subtract are both not infinite or NaN, they have to have the same number of parameters. + * + * Used conventions: + * - \f$ inf - inf = NaN \f$, + * - \f$ -inf - (-inf) = NaN \f$, + * - \f$ NaN - b = NaN \f$, + * - \f$ a - NaN = NaN \f$. + * + * If `std::numeric_limits::has_quiet_NaN` is false, then the returned filtration value will be @ref nan() + * if any operation results in NaN, not only if all operations result in NaN. + * + * @param result First element of the subtraction. + * @param to_subtract Second element of the subtraction. + * @return Entry-wise \f$ result - to\_subtract \f$. + */ + friend One_critical_filtration &operator-=(One_critical_filtration &result, + const One_critical_filtration &to_subtract) { + if (result.empty()) return result; + + if (result.is_nan() || to_subtract.is_nan() || (result.is_plus_inf() && to_subtract.is_plus_inf()) || + (result.is_minus_inf() && to_subtract.is_minus_inf())) { + result = nan(); + return result; + } + if (result.is_plus_inf() || to_subtract.is_minus_inf()) { + result = inf(); + return result; + } + if (result.is_minus_inf() || to_subtract.is_plus_inf()) { + result = minus_inf(); + return result; + } + + GUDHI_CHECK(result.size() == to_subtract.size(), + "Two filtration values with different number of parameters cannot be subtracted."); + + return apply_operation_with_finite_values_(result, to_subtract, subtract_); + } + + /** + * @brief Modifies the first parameters such that an entry at index \f$ i \f$ is equal to + * \f$ result[i] - to\_subtract \f$. + * + * Used conventions: + * - \f$ inf - inf = NaN \f$, + * - \f$ -inf - (-inf) = NaN \f$, + * - \f$ NaN - b = NaN \f$, + * - \f$ a - NaN = NaN \f$. + * + * If `std::numeric_limits::has_quiet_NaN` is false, then the returned filtration value will be @ref nan() + * if any operation results in NaN, not only if all operations result in NaN. + * + * @param result First element of the subtraction. + * @param to_subtract Second element of the subtraction. + * @return Entry-wise \f$ result - to\_subtract \f$. + */ + friend One_critical_filtration &operator-=(One_critical_filtration &result, const T &to_subtract) { + if (result.empty()) return result; + + if (result.is_nan() || is_nan_(to_subtract) || (result.is_plus_inf() && to_subtract == T_inf) || + (result.is_minus_inf() && to_subtract == -T_inf)) { + result = nan(); + return result; + } + if (result.is_plus_inf() || to_subtract == -T_inf) { + result = inf(); + return result; + } + if (result.is_minus_inf() || to_subtract == T_inf) { + result = minus_inf(); + return result; + } + + return apply_scalar_operation_on_finite_value_(result, to_subtract, subtract_); + } + + // addition + /** + * @brief Returns a filtration value such that an entry at index \f$ i \f$ is equal to \f$ a[i] + b[i] \f$. + * If @p a and @p b are both not infinite or NaN, they have to have the same number of parameters. + * + * Used conventions: + * - \f$ NaN + b = NaN \f$, + * - \f$ a + NaN = NaN \f$. + * + * If `std::numeric_limits::has_quiet_NaN` is false, then the returned filtration value will be @ref nan() + * if any operation results in NaN, not only if all operations result in NaN. + * + * @param a First element of the addition. + * @param b Second element of the addition. + * @return Entry-wise \f$ a + b \f$. + */ + friend One_critical_filtration operator+(One_critical_filtration a, const One_critical_filtration &b) { + a += b; + return a; + } + + /** + * @brief Returns a filtration value such that an entry at index \f$ i \f$ is equal to \f$ f[i] + val \f$. + * + * Used conventions: + * - \f$ NaN + b = NaN \f$, + * - \f$ a + NaN = NaN \f$. + * + * If `std::numeric_limits::has_quiet_NaN` is false, then the returned filtration value will be @ref nan() + * if any operation results in NaN, not only if all operations result in NaN. + * + * @param f First element of the addition. + * @param val Second element of the addition. + * @return Entry-wise \f$ f + val \f$. + */ + friend One_critical_filtration operator+(One_critical_filtration f, const T &val) { + f += val; + return f; + } + + /** + * @brief Returns a filtration value such that an entry at index \f$ i \f$ is equal to \f$ val + f[i] \f$. + * + * Used conventions: + * - \f$ NaN + b = NaN \f$, + * - \f$ a + NaN = NaN \f$. + * + * If `std::numeric_limits::has_quiet_NaN` is false, then the returned filtration value will be @ref nan() + * if any operation results in NaN, not only if all operations result in NaN. + * + * @param val First element of the addition. + * @param f Second element of the addition. + * @return Entry-wise \f$ val + f \f$. + */ + friend One_critical_filtration operator+(const T &val, One_critical_filtration f) { + f += val; + return f; + } + + /** + * @brief Modifies the first parameters such that an entry at index \f$ i \f$ is equal to + * \f$ result[i] + to\_add[i] \f$. + * If @p result and @p to_add are both not infinite or NaN, they have to have the same number of parameters. + * + * Used conventions: + * - \f$ NaN + b = NaN \f$, + * - \f$ a + NaN = NaN \f$. + * + * If `std::numeric_limits::has_quiet_NaN` is false, then the returned filtration value will be @ref nan() + * if any operation results in NaN, not only if all operations result in NaN. + * + * @param result First element of the addition. + * @param to_add Second element of the addition. + * @return Entry-wise \f$ result + to\_add \f$. + */ + friend One_critical_filtration &operator+=(One_critical_filtration &result, const One_critical_filtration &to_add) { + if (result.empty()) return result; + + if (result.is_nan() || to_add.is_nan() || (result.is_plus_inf() && to_add.is_minus_inf()) || + (result.is_minus_inf() && to_add.is_plus_inf())) { + result = nan(); + return result; + } + if (result.is_plus_inf() || to_add.is_plus_inf()) { + result = inf(); + return result; + } + if (result.is_minus_inf() || to_add.is_minus_inf()) { + result = minus_inf(); + return result; + } + + GUDHI_CHECK(result.size() == to_add.size(), + "Two filtration values with different number of parameters cannot be added."); + + return apply_operation_with_finite_values_(result, to_add, add_); + } + + /** + * @brief Modifies the first parameters such that an entry at index \f$ i \f$ is equal to + * \f$ result[i] + to\_add \f$. + * + * Used conventions: + * - \f$ NaN + b = NaN \f$, + * - \f$ a + NaN = NaN \f$. + * + * If `std::numeric_limits::has_quiet_NaN` is false, then the returned filtration value will be @ref nan() + * if any operation results in NaN, not only if all operations result in NaN. + * + * @param result First element of the addition. + * @param to_add Second element of the addition. + * @return Entry-wise \f$ result + to\_add \f$. + */ + friend One_critical_filtration &operator+=(One_critical_filtration &result, const T &to_add) { + if (result.empty()) return result; + + if (result.is_nan() || is_nan_(to_add) || (result.is_plus_inf() && to_add == -T_inf) || + (result.is_minus_inf() && to_add == T_inf)) { + result = nan(); + return result; + } + if (result.is_plus_inf() || to_add == T_inf) { + result = inf(); + return result; + } + if (result.is_minus_inf() || to_add == -T_inf) { + result = minus_inf(); + return result; + } + + return apply_scalar_operation_on_finite_value_(result, to_add, add_); + } + + // multiplication + /** + * @brief Returns a filtration value such that an entry at index \f$ i \f$ is equal to \f$ a[i] * b[i] \f$. + * If @p a and @p b are both not infinite or NaN, they have to have the same number of parameters. + * + * Used conventions: + * - \f$ inf * 0 = NaN \f$, + * - \f$ 0 * inf = NaN \f$, + * - \f$ -inf * 0 = NaN \f$, + * - \f$ 0 * -inf = NaN \f$, + * - \f$ NaN * b = NaN \f$, + * - \f$ a * NaN = NaN \f$. + * + * If `std::numeric_limits::has_quiet_NaN` is false, then the returned filtration value will be @ref nan() + * if any operation results in NaN, not only if all operations result in NaN. + * + * @param a First element of the multiplication. + * @param b Second element of the multiplication. + * @return Entry-wise \f$ a * b \f$. + */ + friend One_critical_filtration operator*(One_critical_filtration a, const One_critical_filtration &b) { + a *= b; + return a; + } + + /** + * @brief Returns a filtration value such that an entry at index \f$ i \f$ is equal to \f$ f[i] * val \f$. + * + * Used conventions: + * - \f$ inf * 0 = NaN \f$, + * - \f$ 0 * inf = NaN \f$, + * - \f$ -inf * 0 = NaN \f$, + * - \f$ 0 * -inf = NaN \f$, + * - \f$ NaN * b = NaN \f$, + * - \f$ a * NaN = NaN \f$. + * + * If `std::numeric_limits::has_quiet_NaN` is false, then the returned filtration value will be @ref nan() + * if any operation results in NaN, not only if all operations result in NaN. + * + * @param f First element of the multiplication. + * @param val Second element of the multiplication. + * @return Entry-wise \f$ f * val \f$. + */ + friend One_critical_filtration operator*(One_critical_filtration f, const T &val) { + f *= val; + return f; + } + + /** + * @brief Returns a filtration value such that an entry at index \f$ i \f$ is equal to \f$ val * f[i] \f$. + * + * Used conventions: + * - \f$ inf * 0 = NaN \f$, + * - \f$ 0 * inf = NaN \f$, + * - \f$ -inf * 0 = NaN \f$, + * - \f$ 0 * -inf = NaN \f$, + * - \f$ NaN * b = NaN \f$, + * - \f$ a * NaN = NaN \f$. + * + * If `std::numeric_limits::has_quiet_NaN` is false, then the returned filtration value will be @ref nan() + * if any operation results in NaN, not only if all operations result in NaN. + * + * @param val First element of the multiplication. + * @param f Second element of the multiplication. + * @return Entry-wise \f$ val * f \f$. + */ + friend One_critical_filtration operator*(const T &val, One_critical_filtration f) { + f *= val; + return f; + } + + /** + * @brief Modifies the first parameters such that an entry at index \f$ i \f$ is equal to + * \f$ result[i] * to\_mul[i] \f$. + * If @p result and @p to_mul are both not infinite or NaN, they have to have the same number of parameters. + * + * Used conventions: + * - \f$ inf * 0 = NaN \f$, + * - \f$ 0 * inf = NaN \f$, + * - \f$ -inf * 0 = NaN \f$, + * - \f$ 0 * -inf = NaN \f$, + * - \f$ NaN * b = NaN \f$, + * - \f$ a * NaN = NaN \f$. + * + * If `std::numeric_limits::has_quiet_NaN` is false, then the returned filtration value will be @ref nan() + * if any operation results in NaN, not only if all operations result in NaN. + * + * @param result First element of the multiplication. + * @param to_mul Second element of the multiplication. + * @return Entry-wise \f$ result * to\_mul \f$. + */ + friend One_critical_filtration &operator*=(One_critical_filtration &result, const One_critical_filtration &to_mul) { + if (result.empty()) return result; + + if (result.is_nan() || to_mul.is_nan()) { + result = nan(); + return result; + } + + bool res_is_infinite = result.is_plus_inf() || result.is_minus_inf(); + bool to_mul_is_infinite = to_mul.is_plus_inf() || to_mul.is_minus_inf(); + + if (res_is_infinite && to_mul_is_infinite) { + if (to_mul.is_minus_inf()) { + result[0] = -result[0]; + } + return result; + } + + if (res_is_infinite || to_mul_is_infinite) { + const One_critical_filtration &finite = res_is_infinite ? to_mul : result; + const T infinite = res_is_infinite ? result[0] : to_mul[0]; + result = finite; + return apply_scalar_operation_on_finite_value_(result, infinite, multiply_); + } + + GUDHI_CHECK(result.size() == to_mul.size(), + "Two filtration values with different number of parameters cannot be multiplied."); + + return apply_operation_with_finite_values_(result, to_mul, multiply_); + } + + /** + * @brief Modifies the first parameters such that an entry at index \f$ i \f$ is equal to + * \f$ result[i] * to\_mul \f$. + * + * Used conventions: + * - \f$ inf * 0 = NaN \f$, + * - \f$ 0 * inf = NaN \f$, + * - \f$ -inf * 0 = NaN \f$, + * - \f$ 0 * -inf = NaN \f$, + * - \f$ NaN * b = NaN \f$, + * - \f$ a * NaN = NaN \f$. + * + * If `std::numeric_limits::has_quiet_NaN` is false, then the returned filtration value will be @ref nan() + * if any operation results in NaN, not only if all operations result in NaN. + * + * @param result First element of the multiplication. + * @param to_mul Second element of the multiplication. + * @return Entry-wise \f$ result * to\_mul \f$. + */ + friend One_critical_filtration &operator*=(One_critical_filtration &result, const T &to_mul) { + if (result.empty()) return result; + + if (result.is_nan() || is_nan_(to_mul)) { + result = nan(); + return result; + } + + if (result.is_plus_inf() || result.is_minus_inf()) { + if (to_mul == 0) + result = nan(); + else if (to_mul < 0) + result[0] = -result[0]; + return result; + } + + return apply_scalar_operation_on_finite_value_(result, to_mul, multiply_); + } + + // division + /** + * @brief Returns a filtration value such that an entry at index \f$ i \f$ is equal to \f$ a[i] / b[i] \f$. + * If @p a and @p b are both not infinite or NaN, they have to have the same number of parameters. + * + * Used conventions: + * - \f$ a / 0 = NaN \f$, + * - \f$ inf / inf = NaN \f$, + * - \f$ -inf / inf = NaN \f$, + * - \f$ inf / -inf = NaN \f$, + * - \f$ -inf / -inf = NaN \f$, + * - \f$ NaN / b = NaN \f$, + * - \f$ a / NaN = NaN \f$, + * - \f$ a / inf = 0 \f$, + * - \f$ a / -inf = 0 \f$. + * + * If `std::numeric_limits::has_quiet_NaN` is false, then the returned filtration value will be @ref nan() + * if any operation results in NaN, not only if all operations result in NaN. + * + * @param a First element of the division. + * @param b Second element of the division. + * @return Entry-wise \f$ a / b \f$. + */ + friend One_critical_filtration operator/(One_critical_filtration a, const One_critical_filtration &b) { + a /= b; + return a; + } + + /** + * @brief Returns a filtration value such that an entry at index \f$ i \f$ is equal to \f$ f[i] / val \f$. + * + * Used conventions: + * - \f$ a / 0 = NaN \f$, + * - \f$ inf / inf = NaN \f$, + * - \f$ -inf / inf = NaN \f$, + * - \f$ inf / -inf = NaN \f$, + * - \f$ -inf / -inf = NaN \f$, + * - \f$ NaN / b = NaN \f$, + * - \f$ a / NaN = NaN \f$, + * - \f$ a / inf = 0 \f$, + * - \f$ a / -inf = 0 \f$. + * + * If `std::numeric_limits::has_quiet_NaN` is false, then the returned filtration value will be @ref nan() + * if any operation results in NaN, not only if all operations result in NaN. + * + * @param f First element of the division. + * @param val Second element of the division. + * @return Entry-wise \f$ f / val \f$. + */ + friend One_critical_filtration operator/(One_critical_filtration f, const T &val) { + f /= val; + return f; + } + + /** + * @brief Returns a filtration value such that an entry at index \f$ i \f$ is equal to \f$ val / f[i] \f$. + * + * Used conventions: + * - \f$ a / 0 = NaN \f$, + * - \f$ inf / inf = NaN \f$, + * - \f$ -inf / inf = NaN \f$, + * - \f$ inf / -inf = NaN \f$, + * - \f$ -inf / -inf = NaN \f$, + * - \f$ NaN / b = NaN \f$, + * - \f$ a / NaN = NaN \f$, + * - \f$ a / inf = 0 \f$, + * - \f$ a / -inf = 0 \f$. + * + * If `std::numeric_limits::has_quiet_NaN` is false, then the returned filtration value will be @ref nan() + * if any operation results in NaN, not only if all operations result in NaN. + * + * @param val First element of the division. + * @param f Second element of the division. + * @return Entry-wise \f$ val / f \f$. + */ + friend One_critical_filtration operator/(const T &val, const One_critical_filtration &f) { + if (f.empty()) return f; + if (is_nan_(val) || f.is_nan()) return nan(); + + One_critical_filtration result(f.size(), val); + result /= f; + return result; + } + + /** + * @brief Modifies the first parameters such that an entry at index \f$ i \f$ is equal to + * \f$ result[i] / to\_div[i] \f$. + * If @p result and @p to_div are both not infinite or NaN, they have to have the same number of parameters. + * + * Used conventions: + * - \f$ a / 0 = NaN \f$, + * - \f$ inf / inf = NaN \f$, + * - \f$ -inf / inf = NaN \f$, + * - \f$ inf / -inf = NaN \f$, + * - \f$ -inf / -inf = NaN \f$, + * - \f$ NaN / b = NaN \f$, + * - \f$ a / NaN = NaN \f$, + * - \f$ a / inf = 0 \f$, + * - \f$ a / -inf = 0 \f$. + * + * If `std::numeric_limits::has_quiet_NaN` is false, then the returned filtration value will be @ref nan() + * if any operation results in NaN, not only if all operations result in NaN. + * + * @param result First element of the division. + * @param to_div Second element of the division. + * @return Entry-wise \f$ result / to\_div \f$. + */ + friend One_critical_filtration &operator/=(One_critical_filtration &result, const One_critical_filtration &to_div) { + if (result.empty()) return result; + + bool res_is_infinite = result.is_plus_inf() || result.is_minus_inf(); + bool to_div_is_infinite = to_div.is_plus_inf() || to_div.is_minus_inf(); + + if (result.is_nan() || to_div.is_nan() || (res_is_infinite && to_div_is_infinite)) { + result = nan(); + return result; + } + + if (to_div_is_infinite) { + return apply_scalar_operation_on_finite_value_with_all_nan_possible_(result, to_div[0], divide_); + } + + GUDHI_CHECK(res_is_infinite || result.size() == to_div.size(), + "Two filtration values with different number of parameters cannot be divided."); + + if (res_is_infinite) { + result.resize(to_div.size(), result[0]); + } + + return apply_operation_with_finite_values_(result, to_div, divide_); + } + + /** + * @brief Modifies the first parameters such that an entry at index \f$ i \f$ is equal to + * \f$ result[i] / to\_div \f$. + * + * Used conventions: + * - \f$ a / 0 = NaN \f$, + * - \f$ inf / inf = NaN \f$, + * - \f$ -inf / inf = NaN \f$, + * - \f$ inf / -inf = NaN \f$, + * - \f$ -inf / -inf = NaN \f$, + * - \f$ NaN / b = NaN \f$, + * - \f$ a / NaN = NaN \f$, + * - \f$ a / inf = 0 \f$, + * - \f$ a / -inf = 0 \f$. + * + * If `std::numeric_limits::has_quiet_NaN` is false, then the returned filtration value will be @ref nan() + * if any operation results in NaN, not only if all operations result in NaN. + * + * @param result First element of the division. + * @param to_div Second element of the division. + * @return Entry-wise \f$ result / to\_div \f$. + */ + friend One_critical_filtration &operator/=(One_critical_filtration &result, const T &to_div) { + if (result.empty()) return result; + + bool res_is_infinite = result.is_plus_inf() || result.is_minus_inf(); + bool to_div_is_infinite = to_div == T_inf || to_div == -T_inf; + + if (to_div == 0 || is_nan_(to_div) || result.is_nan() || (res_is_infinite && to_div_is_infinite)) { + result = nan(); + return result; + } + + if (res_is_infinite) { + if (to_div < 0) result[0] = -result[0]; + return result; + } + + return apply_scalar_operation_on_finite_value_with_all_nan_possible_(result, to_div, divide_); + } + + // MODIFIERS + + /** + * @brief Sets the filtration value to the least common upper bound between the current value and the given value. + * + * More formally, it pushes the current filtration value to the cone \f$ \{ y \in \mathbb R^n : y \ge x \} \f$ + * originating in the given filtration value \f$ x \f$. The resulting value corresponds to the intersection of both + * cones: \f$ \mathrm{this} = \min \{ y \in \mathbb R^n : y \ge this \} \cap \{ y \in \mathbb R^n : y \ge x \} \f$. + * + * @param x The target filtration value towards which to push. + */ + void push_to_least_common_upper_bound(const One_critical_filtration &x) { + if (this->is_plus_inf() || this->is_nan() || x.is_nan() || x.is_minus_inf()) return; + if (x.is_plus_inf() || this->is_minus_inf()) { + *this = x; + return; + } + + GUDHI_CHECK(this->num_parameters() == x.num_parameters(), + "A filtration value cannot be pushed to another one with different numbers of parameters."); + + for (std::size_t i = 0; i < x.num_parameters(); i++) + Base::operator[](i) = Base::operator[](i) > x[i] ? Base::operator[](i) : x[i]; + } + + /** + * @brief Sets the filtration value to the greatest common lower bound between the current value and the given value. + * + * More formally, it pulls the current filtration value to the cone \f$ \{ y \in \mathbb R^n : y \le x \} \f$ + * originating in the given filtration value \f$ x \f$. The resulting value corresponds to the intersection of both + * cones: \f$ \mathrm{this} = \min \{ y \in \mathbb R^n : y \le this \} \cap \{ y \in \mathbb R^n : y \le x \} \f$. + * + * @param x The target filtration value towards which to pull. + */ + void pull_to_greatest_common_lower_bound(const One_critical_filtration &x) { + if (x.is_plus_inf() || this->is_nan() || x.is_nan() || this->is_minus_inf()) return; + if (this->is_plus_inf() || x.is_minus_inf()) { + *this = x; + return; + } + + GUDHI_CHECK(this->num_parameters() == x.num_parameters(), + "A filtration value cannot be pulled to another one with different numbers of parameters."); + + for (std::size_t i = 0u; i < x.num_parameters(); i++) + Base::operator[](i) = Base::operator[](i) > x[i] ? x[i] : Base::operator[](i); + } + + /** + * @brief Projects the filtration value into the given grid. If @p coordinate is false, the entries are set to + * the nearest upper bound value with the same parameter in the grid. Otherwise, the entries are set to the indices + * of those nearest upper bound values. + * The grid has to be represented as a vector of ordered ranges of values convertible into `T`. An index + * \f$ i \f$ of the vector corresponds to the same parameter as the index \f$ i \f$ in the filtration value. + * The ranges correspond to the possible values of the parameters, ordered by increasing value, forming therefore + * all together a 2D grid. + * + * @tparam oned_array A range of values convertible into `T` ordered by increasing value. Has to implement + * a begin, end and operator[] method. + * @param grid Vector of @p oned_array with size at least number of filtration parameters. + * @param coordinate If true, the values are set to the coordinates of the projection in the grid. If false, + * the values are set to the values at the coordinates of the projection. + */ + template + void project_onto_grid(const std::vector &grid, bool coordinate = true) { + GUDHI_CHECK(grid.size() >= Base::size(), + "The grid should not be smaller than the number of parameters in the filtration value."); + for (std::size_t parameter = 0u; parameter < Base::size(); ++parameter) { + const auto &filtration = grid[parameter]; + auto d = + std::distance(filtration.begin(), + std::lower_bound(filtration.begin(), filtration.end(), + static_cast(Base::operator[](parameter)))); + Base::operator[](parameter) = coordinate ? static_cast(d) : static_cast(filtration[d]); + } + } + + // FONCTIONNALITIES + + /** + * @brief Computes the scalar product of the given filtration value with the given vector. + * + * @tparam U Arithmetic type of the result. Default value: `T`. + * @param f Filtration value. + * @param x Vector of coefficients. + * @return Scalar product of @p f with @p x. + */ + template + friend U compute_linear_projection(const One_critical_filtration &f, const std::vector &x) { + U projection = 0; + std::size_t size = std::min(x.size(), f.size()); + for (std::size_t i = 0u; i < size; i++) projection += x[i] * static_cast(f[i]); + return projection; + } + + /** + * @brief Computes the norm of the given filtration value. + * + * @param f Filtration value. + * @return The norm of @p f. + */ + friend T compute_norm(const One_critical_filtration &f) { + T out = 0; + for (auto &val : f) out += (val * val); + if constexpr (std::is_integral_v){ + //to avoid Windows issue which don't know how to cast integers for cmath methods + return std::sqrt(static_cast(out)); + } else { + return std::sqrt(out); + } + } + + /** + * @brief Computes the euclidean distance from the first parameter to the second parameter. + * + * @param f Start filtration value. + * @param other End filtration value. + * @return Euclidean distance between @p f and @p other. + */ + friend T compute_euclidean_distance_to(const One_critical_filtration &f, const One_critical_filtration &other) { + T out = 0; + for (std::size_t i = 0u; i < other.size(); i++) { + out += (f[i] - other[i]) * (f[i] - other[i]); + } + if constexpr (std::is_integral_v){ + //to avoid Windows issue which don't know how to cast integers for cmath methods + return std::sqrt(static_cast(out)); + } else { + return std::sqrt(out); + } + } + + /** + * @brief Computes the coordinates in the given grid, corresponding to the nearest upper bounds of the entries + * in the given filtration value. + * The grid has to be represented as a vector of vectors of ordered values convertible into `out_type`. An index + * \f$ i \f$ of the vector corresponds to the same parameter as the index \f$ i \f$ in the filtration value. + * The inner vectors correspond to the possible values of the parameters, ordered by increasing value, + * forming therefore all together a 2D grid. + * + * @tparam out_type Signed arithmetic type. Default value: std::int32_t. + * @tparam U Type which is convertible into `out_type`. + * @param f Filtration value to project. + * @param grid Vector of vectors to project into. + * @return Filtration value \f$ out \f$ whose entry correspond to the indices of the projected values. That is, + * the projection of \f$ f[i] \f$ is \f$ grid[i][out[i]] \f$. + */ + template + friend One_critical_filtration compute_coordinates_in_grid(One_critical_filtration f, + const std::vector > &grid) { + // TODO: by replicating the code of "project_onto_grid", this could be done with just one copy + // instead of two. But it is not clear if it is really worth it, i.e., how much the change in type is really + // necessary in the use cases. To see later. + f.project_onto_grid(grid); + if constexpr (std::is_same_v) { + return f; + } else { + return f.as_type(); + } + } + + /** + * @brief Computes the values in the given grid corresponding to the coordinates given by the given filtration + * value. That is, if \f$ out \f$ is the result, \f$ out[i] = grid[i][f[i]] \f$. Assumes therefore, that the + * values stored in the filtration value corresponds to indices existing in the given grid. + * + * @tparam U Signed arithmetic type. + * @param f Filtration value storing coordinates compatible with `grid`. + * @param grid Vector of vector. + * @return Filtration value \f$ out \f$ whose entry correspond to \f$ out[i] = grid[i][f[i]] \f$. + */ + template + friend One_critical_filtration evaluate_coordinates_in_grid(const One_critical_filtration &f, + const std::vector > &grid) { + One_critical_filtration pushed_value(f.size()); + + GUDHI_CHECK(grid.size() == f.size(), + "The size of the grid should correspond to the number of parameters in the filtration value."); + + U grid_inf = One_critical_filtration::T_inf; + + for (std::size_t parameter = 0u; parameter < grid.size(); ++parameter) { + const auto &filtration = grid[parameter]; + const auto &c = f[parameter]; + pushed_value[parameter] = c == f.T_inf ? grid_inf : filtration[c]; + } + return pushed_value; + } + + // UTILITIES + + /** + * @brief Outstream operator. + */ + friend std::ostream &operator<<(std::ostream &stream, const One_critical_filtration &f) { + if (f.is_plus_inf()) { + stream << "[inf, ..., inf]"; + return stream; + } + if (f.is_minus_inf()) { + stream << "[-inf, ..., -inf]"; + return stream; + } + if (f.is_nan()) { + stream << "[NaN]"; + return stream; + } + if (f.empty()) { + stream << "[]"; + return stream; + } + stream << "["; + for (std::size_t i = 0; i < f.size() - 1; i++) { + stream << f[i] << ", "; + } + stream << f.back(); + stream << "]"; + return stream; + } + + public: + /** + * @brief Infinity value of an entry of the filtration value. + */ + constexpr static const T T_inf = + std::numeric_limits::has_infinity ? std::numeric_limits::infinity() : std::numeric_limits::max(); + + /** + * @brief Indicates if the class manages multi-critical filtration values. + */ + constexpr static bool is_multi_critical = false; + + private: + static bool is_nan_(T val){ + if constexpr (std::is_integral_v){ + //to avoid Windows issue which don't know how to cast integers for cmath methods + return false; + } else { + return std::isnan(val); + } + } + + constexpr static bool subtract_(T &v1, T v2) { return add_(v1, -v2); } + + constexpr static bool add_(T &v1, T v2) { + if (is_nan_(v1) || is_nan_(v2) || (v1 == T_inf && v2 == -T_inf) || (v1 == -T_inf && v2 == T_inf)) { + v1 = std::numeric_limits::quiet_NaN(); + return false; + } + if (v1 == T_inf || v1 == -T_inf) { + return true; + } + if (v2 == T_inf || v2 == -T_inf) { + v1 = v2; + return true; + } + + v1 += v2; + return true; + } + + constexpr static bool multiply_(T &v1, T v2) { + bool v1_is_infinite = v1 == T_inf || v1 == -T_inf; + bool v2_is_infinite = v2 == T_inf || v2 == -T_inf; + + if (is_nan_(v1) || is_nan_(v2) || (v1_is_infinite && v2 == 0) || (v1 == 0 && v2_is_infinite)) { + v1 = std::numeric_limits::quiet_NaN(); + return false; + } + + if ((v1 == T_inf && v2 > 0) || (v1 == -T_inf && v2 < 0) || (v1 < 0 && v2 == -T_inf) || (v1 > 0 && v2 == T_inf)) { + v1 = T_inf; + return true; + } + + if ((v1 == T_inf && v2 < 0) || (v1 == -T_inf && v2 > 0) || (v1 > 0 && v2 == -T_inf) || (v1 < 0 && v2 == T_inf)) { + v1 = -T_inf; + return true; + } + + v1 *= v2; + return true; + } + + constexpr static bool divide_(T &v1, T v2) { + bool v1_is_infinite = v1 == T_inf || v1 == -T_inf; + bool v2_is_infinite = v2 == T_inf || v2 == -T_inf; + + if (is_nan_(v1) || is_nan_(v2) || v2 == 0 || (v1_is_infinite && v2_is_infinite)) { + v1 = std::numeric_limits::quiet_NaN(); + return false; + } + + if (v1 == 0 || (v1_is_infinite && v2 > 0)) return true; + + if (v1_is_infinite && v2 < 0) { + v1 = -v1; + return true; + } + + if (v2_is_infinite) { + v1 = 0; + return true; + } + + v1 /= v2; + return true; + } + + constexpr static bool update_sign_(T toComp, int &sign) { + if (toComp == T_inf) { + if (sign == 0) + sign = 1; + else if (sign == -1) + return false; + } else if (toComp == -T_inf) { + if (sign == 0) + sign = -1; + else if (sign == 1) + return false; + } else { + return false; + } + + return true; + } + + template + static One_critical_filtration &apply_operation_with_finite_values_(One_critical_filtration &result, + const One_critical_filtration &to_operate, + F &&operate) { + bool allSameInf = true; + bool allNaN = true; + int sign = 0; + for (auto i = 0u; i < result.size(); ++i) { + if (operate(result[i], to_operate[i])) { + allNaN = false; + } else { + if constexpr (!std::numeric_limits::has_quiet_NaN) { + result = nan(); + return result; + } + } + if (allSameInf) allSameInf = update_sign_(result[i], sign); + } + + if (allSameInf) result = (sign == 1 ? inf() : minus_inf()); + if (allNaN) result = nan(); + + return result; + } + + template + static One_critical_filtration &apply_scalar_operation_on_finite_value_(One_critical_filtration &result, + const T &to_operate, F &&operate) { + for (auto &val : result) { + if constexpr (std::numeric_limits::has_quiet_NaN) { + operate(val, to_operate); + } else { + if (!operate(val, to_operate)) { + result = nan(); + return result; + } + } + } + + return result; + } + + template + static One_critical_filtration &apply_scalar_operation_on_finite_value_with_all_nan_possible_( + One_critical_filtration &result, const T &to_operate, F &&operate) { + bool allNaN = true; + + for (auto &val : result) { + if (operate(val, to_operate)) { + allNaN = false; + } else { + if constexpr (!std::numeric_limits::has_quiet_NaN) { + result = nan(); + return result; + } + } + } + if (allNaN) result = nan(); + + return result; + } +}; + +} // namespace Gudhi::multi_filtration + +namespace std { + +template +class numeric_limits > { + public: + static constexpr bool has_infinity = true; + + static constexpr Gudhi::multi_filtration::One_critical_filtration infinity() noexcept { + return Gudhi::multi_filtration::One_critical_filtration::inf(); + }; + + // non-standard + static constexpr Gudhi::multi_filtration::One_critical_filtration minus_infinity() noexcept { + return Gudhi::multi_filtration::One_critical_filtration::minus_inf(); + }; + + static constexpr Gudhi::multi_filtration::One_critical_filtration max() noexcept(false) { + throw std::logic_error( + "The maximal value cannot be represented with no finite numbers of parameters." + "Use `max(number_of_parameters)` instead"); + }; + + // non-standard, so I don't want to define a default value. + static constexpr Gudhi::multi_filtration::One_critical_filtration max(unsigned int n) noexcept { + return Gudhi::multi_filtration::One_critical_filtration(n, std::numeric_limits::max()); + }; + + static constexpr Gudhi::multi_filtration::One_critical_filtration quiet_NaN() noexcept { + return Gudhi::multi_filtration::One_critical_filtration::nan(); + }; +}; + +} // namespace std + +#endif // ONE_CRITICAL_FILTRATIONS_H_ diff --git a/src/Multi_filtration/test/CMakeLists.txt b/src/Multi_filtration/test/CMakeLists.txt new file mode 100644 index 0000000000..2167fa080b --- /dev/null +++ b/src/Multi_filtration/test/CMakeLists.txt @@ -0,0 +1,7 @@ +include(GUDHI_boost_test) + +add_executable_with_targets(Multi_filtration_onecritical_unit_test multifiltration_onecritical_unit_test.cpp TBB::tbb) +gudhi_add_boost_test(Multi_filtration_onecritical_unit_test) + +add_executable_with_targets(Multi_filtration_multicritical_unit_test multifiltration_multicritical_unit_test.cpp TBB::tbb) +gudhi_add_boost_test(Multi_filtration_multicritical_unit_test) diff --git a/src/Multi_filtration/test/multifiltration_multicritical_unit_test.cpp b/src/Multi_filtration/test/multifiltration_multicritical_unit_test.cpp new file mode 100644 index 0000000000..5c18498161 --- /dev/null +++ b/src/Multi_filtration/test/multifiltration_multicritical_unit_test.cpp @@ -0,0 +1,473 @@ +/* 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): Hannah Schreiber + * + * Copyright (C) 2024 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#include +#include +#include + +#define BOOST_TEST_DYN_LINK +#define BOOST_TEST_MODULE "multi_filtration" +#include +#include + +#include + +using Gudhi::multi_filtration::Multi_critical_filtration; + +typedef boost::mpl::list list_of_tested_variants; + +BOOST_AUTO_TEST_CASE_TEMPLATE(multi_critical_filtration_constructors, T, list_of_tested_variants) +{ + Multi_critical_filtration f01; + BOOST_CHECK(f01.num_parameters() == 1); + BOOST_CHECK(f01.num_generators() == 1); + BOOST_CHECK(f01[0][0] == -Multi_critical_filtration::Generator::T_inf); + Multi_critical_filtration f02; + BOOST_CHECK(f02.num_parameters() == 1); + BOOST_CHECK(f02.num_generators() == 1); + BOOST_CHECK(f02[0][0] == Multi_critical_filtration::Generator::T_inf); + + Multi_critical_filtration f1(3); + BOOST_CHECK(f1.num_parameters() == 3); + BOOST_CHECK(f1.num_generators() == 1); + BOOST_CHECK(f1[0][0] == -Multi_critical_filtration::Generator::T_inf); + BOOST_CHECK(f1[0][1] == -Multi_critical_filtration::Generator::T_inf); + BOOST_CHECK(f1[0][2] == -Multi_critical_filtration::Generator::T_inf); + + Multi_critical_filtration f2(3, 0); + BOOST_CHECK(f2.num_parameters() == 3); + BOOST_CHECK(f2.num_generators() == 1); + BOOST_CHECK(f2[0][0] == 0); + BOOST_CHECK(f2[0][1] == 0); + BOOST_CHECK(f2[0][2] == 0); + + Multi_critical_filtration f3({0, 1, 2}); + BOOST_CHECK(f3.num_parameters() == 3); + BOOST_CHECK(f3.num_generators() == 1); + BOOST_CHECK(f3[0][0] == 0); + BOOST_CHECK(f3[0][1] == 1); + BOOST_CHECK(f3[0][2] == 2); + + std::vector v{0, 1, 2}; + Multi_critical_filtration f41(v); + BOOST_CHECK(f41.num_parameters() == 3); + BOOST_CHECK(f41.num_generators() == 1); + BOOST_CHECK(f41[0][0] == 0); + BOOST_CHECK(f41[0][1] == 1); + BOOST_CHECK(f41[0][2] == 2); + + Multi_critical_filtration f5(v.begin(), v.end()); + BOOST_CHECK(f5.num_generators() == 1); + BOOST_CHECK(f5.num_parameters() == 3); + BOOST_CHECK(f5[0][0] == 0); + BOOST_CHECK(f5[0][1] == 1); + BOOST_CHECK(f5[0][2] == 2); + + Multi_critical_filtration f42(std::move(v)); + BOOST_CHECK(f42.num_parameters() == 3); + BOOST_CHECK(f42.num_generators() == 1); + BOOST_CHECK(f42[0][0] == 0); + BOOST_CHECK(f42[0][1] == 1); + BOOST_CHECK(f42[0][2] == 2); + + std::vector::Generator> v2{{0, 1, 2}, {3, 4, 5}}; + Multi_critical_filtration f6(std::move(v2)); + BOOST_CHECK(f6.num_parameters() == 3); + BOOST_CHECK(f6.num_generators() == 2); + BOOST_CHECK(f6[0][0] == 0); + BOOST_CHECK(f6[0][1] == 1); + BOOST_CHECK(f6[0][2] == 2); + BOOST_CHECK(f6[1][0] == 3); + BOOST_CHECK(f6[1][1] == 4); + BOOST_CHECK(f6[1][2] == 5); + + Multi_critical_filtration f7(f6); + BOOST_CHECK(f7.num_parameters() == 3); + BOOST_CHECK(f7.num_generators() == 2); + BOOST_CHECK(f7[0][0] == 0); + BOOST_CHECK(f7[0][1] == 1); + BOOST_CHECK(f7[0][2] == 2); + BOOST_CHECK(f7[1][0] == 3); + BOOST_CHECK(f7[1][1] == 4); + BOOST_CHECK(f7[1][2] == 5); + + f01 = f5; + BOOST_CHECK(f01.num_generators() == 1); + BOOST_CHECK(f01.num_parameters() == 3); + BOOST_CHECK(f01[0][0] == 0); + BOOST_CHECK(f01[0][1] == 1); + BOOST_CHECK(f01[0][2] == 2); +} + +BOOST_AUTO_TEST_CASE_TEMPLATE(multi_critical_filtration_utilities, T, list_of_tested_variants) +{ + Multi_critical_filtration f({0, 1, 2}); + bool test = std::is_same_v; + BOOST_CHECK(test); + + Multi_critical_filtration f2 = f.template as_type(); + test = std::is_same_v; + BOOST_CHECK(test); + BOOST_CHECK(f2.num_generators() == 1); + BOOST_CHECK(f2.num_parameters() == 3); + BOOST_CHECK(f2[0][0] == 0.); + BOOST_CHECK(f2[0][1] == 1.); + BOOST_CHECK(f2[0][2] == 2.); + + BOOST_CHECK(!f.is_plus_inf()); + BOOST_CHECK(!f.is_minus_inf()); + BOOST_CHECK(!f.is_nan()); + BOOST_CHECK(f.is_finite()); + + Multi_critical_filtration f31; + BOOST_CHECK(!f31.is_plus_inf()); + BOOST_CHECK(f31.is_minus_inf()); + BOOST_CHECK(!f31.is_nan()); + BOOST_CHECK(!f31.is_finite()); + + Multi_critical_filtration f32; + BOOST_CHECK(f32.is_plus_inf()); + BOOST_CHECK(!f32.is_minus_inf()); + BOOST_CHECK(!f32.is_nan()); + BOOST_CHECK(!f32.is_finite()); + + //{-inf, -inf, -inf} is considered finite as the user is supposed to updates the values to something else + //the idea is just to reserve space and to give the possibility to use `f4[i] =` + //if the value should really be -inf, use `f4(1)` or `f4 = minus_inf()` instead. + Multi_critical_filtration f4(3); + BOOST_CHECK(!f4.is_plus_inf()); + BOOST_CHECK(!f4.is_minus_inf()); + BOOST_CHECK(!f4.is_nan()); + BOOST_CHECK(f4.is_finite()); + + Multi_critical_filtration f5(1); + BOOST_CHECK(!f5.is_plus_inf()); + BOOST_CHECK(f5.is_minus_inf()); + BOOST_CHECK(!f5.is_nan()); + BOOST_CHECK(!f5.is_finite()); + + auto f6 = Multi_critical_filtration::nan(); + BOOST_CHECK(!f6.is_plus_inf()); + BOOST_CHECK(!f6.is_minus_inf()); + BOOST_CHECK(f6.is_nan()); + BOOST_CHECK(!f6.is_finite()); + + auto f7 = Multi_critical_filtration::inf(); + BOOST_CHECK(f7.is_plus_inf()); + BOOST_CHECK(!f7.is_minus_inf()); + BOOST_CHECK(!f7.is_nan()); + BOOST_CHECK(!f7.is_finite()); +} + +BOOST_AUTO_TEST_CASE_TEMPLATE(multi_critical_filtration_comparators, T, list_of_tested_variants) +{ + Multi_critical_filtration f1({{-1, 4, 5}, {0, 1, 2}}); + Multi_critical_filtration f2({{-5, 0, 0}, {-2, -5, -1}, {-2, 0, 1}}); + Multi_critical_filtration f3({4, 5, 6}); + Multi_critical_filtration f4({{-4, 5, 6}, {0, 0, 1}}); + + BOOST_CHECK(!(f1 < f1)); + BOOST_CHECK(!(f1 < f2)); + BOOST_CHECK(f1 < f3); + BOOST_CHECK(!(f1 < f4)); + BOOST_CHECK(!(f1 < Multi_critical_filtration::nan())); + BOOST_CHECK(f1 < Multi_critical_filtration::inf()); + BOOST_CHECK(!(f1 < Multi_critical_filtration::minus_inf())); + + BOOST_CHECK(f1 <= f1); + BOOST_CHECK(!(f1 <= f2)); + BOOST_CHECK(f1 <= f3); + BOOST_CHECK(!(f1 <= f4)); + BOOST_CHECK(!(f1 <= Multi_critical_filtration::nan())); + BOOST_CHECK(f1 <= Multi_critical_filtration::inf()); + BOOST_CHECK(!(f1 <= Multi_critical_filtration::minus_inf())); + + BOOST_CHECK(!(f1 > f1)); + BOOST_CHECK(f1 > f2); + BOOST_CHECK(!(f1 > f3)); + BOOST_CHECK(!(f1 > f4)); + BOOST_CHECK(!(f1 > Multi_critical_filtration::nan())); + BOOST_CHECK(!(f1 > Multi_critical_filtration::inf())); + BOOST_CHECK(f1 > Multi_critical_filtration::minus_inf()); + + BOOST_CHECK(f1 >= f1); + BOOST_CHECK(f1 >= f2); + BOOST_CHECK(!(f1 >= f3)); + BOOST_CHECK(!(f1 >= f4)); + BOOST_CHECK(!(f1 >= Multi_critical_filtration::nan())); + BOOST_CHECK(!(f1 >= Multi_critical_filtration::inf())); + BOOST_CHECK(f1 >= Multi_critical_filtration::minus_inf()); + + BOOST_CHECK(f1 == f1); + BOOST_CHECK(!(f1 == f2)); + BOOST_CHECK(!(f1 == f3)); + BOOST_CHECK(!(f1 == f4)); + BOOST_CHECK(!(f1 == Multi_critical_filtration::nan())); + BOOST_CHECK(!(f1 == Multi_critical_filtration::inf())); + BOOST_CHECK(!(f1 == Multi_critical_filtration::minus_inf())); + + BOOST_CHECK(!(f1 != f1)); + BOOST_CHECK(f1 != f2); + BOOST_CHECK(f1 != f3); + BOOST_CHECK(f1 != f4); + BOOST_CHECK(f1 != Multi_critical_filtration::nan()); + BOOST_CHECK(f1 != Multi_critical_filtration::inf()); + BOOST_CHECK(f1 != Multi_critical_filtration::minus_inf()); +} + +BOOST_AUTO_TEST_CASE_TEMPLATE(multi_critical_filtration_modifiers, T, list_of_tested_variants) +{ + Multi_critical_filtration f({{-3, 1, 7}, {0, 1, 2}}); + BOOST_CHECK_EQUAL(f[0][0], -3); + BOOST_CHECK_EQUAL(f[0][1], 1); + BOOST_CHECK_EQUAL(f[0][2], 7); + BOOST_CHECK_EQUAL(f[1][0], 0); + BOOST_CHECK_EQUAL(f[1][1], 1); + BOOST_CHECK_EQUAL(f[1][2], 2); + + f.push_to_least_common_upper_bound({-1, 5, 6}); + BOOST_CHECK_EQUAL(f[0][0], -1); + BOOST_CHECK_EQUAL(f[0][1], 5); + BOOST_CHECK_EQUAL(f[0][2], 7); + BOOST_CHECK_EQUAL(f[1][0], 0); + BOOST_CHECK_EQUAL(f[1][1], 5); + BOOST_CHECK_EQUAL(f[1][2], 6); + + f.push_to_least_common_upper_bound({-1, -5, -6}); + BOOST_CHECK_EQUAL(f[0][0], -1); + BOOST_CHECK_EQUAL(f[0][1], 5); + BOOST_CHECK_EQUAL(f[0][2], 7); + BOOST_CHECK_EQUAL(f[1][0], 0); + BOOST_CHECK_EQUAL(f[1][1], 5); + BOOST_CHECK_EQUAL(f[1][2], 6); + + f.push_to_least_common_upper_bound(Multi_critical_filtration::Generator::minus_inf()); + BOOST_CHECK_EQUAL(f[0][0], -1); + BOOST_CHECK_EQUAL(f[0][1], 5); + BOOST_CHECK_EQUAL(f[0][2], 7); + BOOST_CHECK_EQUAL(f[1][0], 0); + BOOST_CHECK_EQUAL(f[1][1], 5); + BOOST_CHECK_EQUAL(f[1][2], 6); + + f.push_to_least_common_upper_bound(Multi_critical_filtration::Generator::inf()); + BOOST_CHECK(f.is_plus_inf()); + + f.push_to_least_common_upper_bound(Multi_critical_filtration::Generator::nan()); + BOOST_CHECK(f.is_plus_inf()); + + f.pull_to_greatest_common_lower_bound({-1, 5, 6}); + BOOST_CHECK_EQUAL(f[0][0], -1); + BOOST_CHECK_EQUAL(f[0][1], 5); + BOOST_CHECK_EQUAL(f[0][2], 6); + + f.pull_to_greatest_common_lower_bound({1, 8, 9}); + BOOST_CHECK_EQUAL(f[0][0], -1); + BOOST_CHECK_EQUAL(f[0][1], 5); + BOOST_CHECK_EQUAL(f[0][2], 6); + + f.pull_to_greatest_common_lower_bound(Multi_critical_filtration::Generator::inf()); + BOOST_CHECK_EQUAL(f[0][0], -1); + BOOST_CHECK_EQUAL(f[0][1], 5); + BOOST_CHECK_EQUAL(f[0][2], 6); + + f.pull_to_greatest_common_lower_bound(Multi_critical_filtration::Generator::minus_inf()); + BOOST_CHECK(f.is_minus_inf()); + + f.pull_to_greatest_common_lower_bound(Multi_critical_filtration::Generator::nan()); + BOOST_CHECK(f.is_minus_inf()); + + std::vector > grid = {{0, 2, 4, 8}, {0, 3, 6, 9}, {0, 4, 8, 16}}; + + f.push_to_least_common_upper_bound({1, 7, 5}); + f.project_onto_grid(grid, true); + BOOST_CHECK_EQUAL(f[0][0], 1); + BOOST_CHECK_EQUAL(f[0][1], 3); + BOOST_CHECK_EQUAL(f[0][2], 2); + + f.push_to_least_common_upper_bound({1, 7, 5}); + f.project_onto_grid(grid, false); + BOOST_CHECK_EQUAL(f[0][0], 2); + BOOST_CHECK_EQUAL(f[0][1], 9); + BOOST_CHECK_EQUAL(f[0][2], 8); +} + +BOOST_AUTO_TEST_CASE_TEMPLATE(multi_critical_filtration_add, T, list_of_tested_variants) +{ + Multi_critical_filtration f({0, 1, 2}); + BOOST_CHECK_EQUAL(f.num_generators(), 1); + BOOST_CHECK_EQUAL(f.num_parameters(), 3); + BOOST_CHECK_EQUAL(f[0][0], 0); + BOOST_CHECK_EQUAL(f[0][1], 1); + BOOST_CHECK_EQUAL(f[0][2], 2); + + bool res = f.add_generator({-3, 1, 7}); + BOOST_CHECK(res); + BOOST_CHECK_EQUAL(f.num_generators(), 2); + BOOST_CHECK_EQUAL(f.num_parameters(), 3); + BOOST_CHECK_EQUAL(f[0][0], -3); + BOOST_CHECK_EQUAL(f[0][1], 1); + BOOST_CHECK_EQUAL(f[0][2], 7); + BOOST_CHECK_EQUAL(f[1][0], 0); + BOOST_CHECK_EQUAL(f[1][1], 1); + BOOST_CHECK_EQUAL(f[1][2], 2); + + res = f.add_generator({-1, -2, -3}); + BOOST_CHECK(res); + BOOST_CHECK_EQUAL(f.num_generators(), 2); + BOOST_CHECK_EQUAL(f.num_parameters(), 3); + BOOST_CHECK_EQUAL(f[0][0], -3); + BOOST_CHECK_EQUAL(f[0][1], 1); + BOOST_CHECK_EQUAL(f[0][2], 7); + BOOST_CHECK_EQUAL(f[1][0], -1); + BOOST_CHECK_EQUAL(f[1][1], -2); + BOOST_CHECK_EQUAL(f[1][2], -3); + + res = f.add_generator({8, 9, 10}); + BOOST_CHECK(!res); + BOOST_CHECK_EQUAL(f.num_generators(), 2); + BOOST_CHECK_EQUAL(f.num_parameters(), 3); + BOOST_CHECK_EQUAL(f[0][0], -3); + BOOST_CHECK_EQUAL(f[0][1], 1); + BOOST_CHECK_EQUAL(f[0][2], 7); + BOOST_CHECK_EQUAL(f[1][0], -1); + BOOST_CHECK_EQUAL(f[1][1], -2); + BOOST_CHECK_EQUAL(f[1][2], -3); + + res = f.add_generator(Multi_critical_filtration::Generator::inf()); + BOOST_CHECK(!res); + BOOST_CHECK_EQUAL(f.num_generators(), 2); + BOOST_CHECK_EQUAL(f.num_parameters(), 3); + BOOST_CHECK_EQUAL(f[0][0], -3); + BOOST_CHECK_EQUAL(f[0][1], 1); + BOOST_CHECK_EQUAL(f[0][2], 7); + BOOST_CHECK_EQUAL(f[1][0], -1); + BOOST_CHECK_EQUAL(f[1][1], -2); + BOOST_CHECK_EQUAL(f[1][2], -3); + + res = f.add_generator(Multi_critical_filtration::Generator::nan()); + BOOST_CHECK(!res); + BOOST_CHECK_EQUAL(f.num_generators(), 2); + BOOST_CHECK_EQUAL(f.num_parameters(), 3); + BOOST_CHECK_EQUAL(f[0][0], -3); + BOOST_CHECK_EQUAL(f[0][1], 1); + BOOST_CHECK_EQUAL(f[0][2], 7); + BOOST_CHECK_EQUAL(f[1][0], -1); + BOOST_CHECK_EQUAL(f[1][1], -2); + BOOST_CHECK_EQUAL(f[1][2], -3); + + res = f.add_generator(Multi_critical_filtration::Generator::minus_inf()); + BOOST_CHECK(res); + BOOST_CHECK_EQUAL(f.num_generators(), 1); + BOOST_CHECK_EQUAL(f.num_parameters(), 1); + BOOST_CHECK_EQUAL(f[0][0], -Multi_critical_filtration::Generator::T_inf); + + std::vector::Generator> v{ + {0, 1, 2}, + typename Multi_critical_filtration::Generator(0), + Multi_critical_filtration::Generator::inf(), + {0, 1, 2}, + Multi_critical_filtration::Generator::nan(), + typename Multi_critical_filtration::Generator(0), + Multi_critical_filtration::Generator::minus_inf()}; + + Multi_critical_filtration f2(v); + f2.remove_empty_generators(false); + BOOST_CHECK_EQUAL(f2[0][0], 0); + BOOST_CHECK_EQUAL(f2[0][1], 1); + BOOST_CHECK_EQUAL(f2[0][2], 2); + BOOST_CHECK_EQUAL(f2[1][0], 0); + BOOST_CHECK_EQUAL(f2[1][1], 1); + BOOST_CHECK_EQUAL(f2[1][2], 2); + BOOST_CHECK(f2[2].is_minus_inf()); + BOOST_CHECK(f2[3].is_plus_inf()); + if constexpr (std::numeric_limits::has_quiet_NaN){ + BOOST_CHECK_EQUAL(f2.num_generators(), 5); + BOOST_CHECK(f2[4].is_nan()); + } else { + BOOST_CHECK_EQUAL(f2.num_generators(), 4); + } + + Multi_critical_filtration f3(v); + f3.remove_empty_generators(true); + BOOST_CHECK_EQUAL(f3[0][0], 0); + BOOST_CHECK_EQUAL(f3[0][1], 1); + BOOST_CHECK_EQUAL(f3[0][2], 2); + BOOST_CHECK_EQUAL(f3[1][0], 0); + BOOST_CHECK_EQUAL(f3[1][1], 1); + BOOST_CHECK_EQUAL(f3[1][2], 2); + if constexpr (std::numeric_limits::has_quiet_NaN){ + BOOST_CHECK_EQUAL(f3.num_generators(), 3); + BOOST_CHECK(f3[2].is_nan()); + } else { + BOOST_CHECK_EQUAL(f3.num_generators(), 2); + } + + Multi_critical_filtration f4(v); + f4.simplify(); + BOOST_CHECK_EQUAL(f4.num_generators(), 1); + BOOST_CHECK(f4[0].is_minus_inf()); + + v.pop_back(); + Multi_critical_filtration f5(v); + f5.simplify(); + BOOST_CHECK_EQUAL(f5.num_generators(), 1); + BOOST_CHECK_EQUAL(f5[0][0], 0); + BOOST_CHECK_EQUAL(f5[0][1], 1); + BOOST_CHECK_EQUAL(f5[0][2], 2); + + Multi_critical_filtration f6 = Multi_critical_filtration::minus_inf(); + bool change = f6.add_generator(Multi_critical_filtration::inf()); + BOOST_CHECK(change); + BOOST_CHECK(f6.is_plus_inf()); +} + +BOOST_AUTO_TEST_CASE_TEMPLATE(multi_critical_filtration_friends, T, list_of_tested_variants) +{ + Multi_critical_filtration f({{0, 1, 2}, {2, 0, 4}}); + + BOOST_CHECK_EQUAL(compute_linear_projection(f, {2,3,5,9}), 13); + BOOST_CHECK(factorize_below(f) == typename Multi_critical_filtration::Generator({0, 0, 2})); + BOOST_CHECK(factorize_above(f) == typename Multi_critical_filtration::Generator({2, 1, 4})); + + Multi_critical_filtration f2({{0, 1, 2}, {2, 0, 4}}); + BOOST_CHECK_EQUAL(compute_linear_projection(f2, {2,3,5,9}), 24); + + f[0] = {1, 7, 5}; + + std::vector > grid = {{0, 2, 4, 8}, {0, 3, 6, 9}, {0, 4, 8, 16}}; + auto res = compute_coordinates_in_grid(f, grid); + BOOST_CHECK_EQUAL(res[0][0], 1); + BOOST_CHECK_EQUAL(res[0][1], 3); + BOOST_CHECK_EQUAL(res[0][2], 2); + BOOST_CHECK_EQUAL(res[1][0], 1); + BOOST_CHECK_EQUAL(res[1][1], 0); + BOOST_CHECK_EQUAL(res[1][2], 1); + + res = evaluate_coordinates_in_grid(res, grid); + BOOST_CHECK_EQUAL(res[0][0], 2); + BOOST_CHECK_EQUAL(res[0][1], 0); + BOOST_CHECK_EQUAL(res[0][2], 4); +} + +BOOST_AUTO_TEST_CASE_TEMPLATE(multi_critical_filtration_numerical_limits, T, list_of_tested_variants) +{ + BOOST_CHECK(std::numeric_limits >::has_infinity); + BOOST_CHECK(std::numeric_limits >::infinity().is_plus_inf()); + BOOST_CHECK(std::numeric_limits >::quiet_NaN().is_nan()); + BOOST_CHECK_THROW(std::numeric_limits >::max(), std::logic_error); + auto max = std::numeric_limits >::max(2, 3); + BOOST_CHECK_EQUAL(max[0][0], std::numeric_limits::max()); + BOOST_CHECK_EQUAL(max[0][1], std::numeric_limits::max()); + BOOST_CHECK_EQUAL(max[0][2], std::numeric_limits::max()); + BOOST_CHECK_EQUAL(max[1][0], std::numeric_limits::max()); + BOOST_CHECK_EQUAL(max[1][1], std::numeric_limits::max()); + BOOST_CHECK_EQUAL(max[1][2], std::numeric_limits::max()); +} + diff --git a/src/Multi_filtration/test/multifiltration_onecritical_unit_test.cpp b/src/Multi_filtration/test/multifiltration_onecritical_unit_test.cpp new file mode 100644 index 0000000000..b7924405be --- /dev/null +++ b/src/Multi_filtration/test/multifiltration_onecritical_unit_test.cpp @@ -0,0 +1,504 @@ +/* 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): Hannah Schreiber + * + * Copyright (C) 2024 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#include +#include + +#define BOOST_TEST_DYN_LINK +#define BOOST_TEST_MODULE "multi_filtration" +#include +#include + +#include + +using Gudhi::multi_filtration::One_critical_filtration; + +typedef boost::mpl::list list_of_tested_variants; + +BOOST_AUTO_TEST_CASE_TEMPLATE(one_critical_filtration_constructors, T, list_of_tested_variants) +{ + One_critical_filtration f; + BOOST_CHECK(!f.empty()); + BOOST_CHECK(f.num_parameters() == 1); + + One_critical_filtration f1(3); + BOOST_CHECK(f1.size() == 3); + BOOST_CHECK(f1.num_parameters() == 3); + BOOST_CHECK(f1[0] == -One_critical_filtration::T_inf); + + One_critical_filtration f2(3, 0); + BOOST_CHECK(f2.size() == 3); + BOOST_CHECK(f2.num_parameters() == 3); + BOOST_CHECK(f2[0] == 0); + + One_critical_filtration f3({0, 1, 2}); + BOOST_CHECK(f3.size() == 3); + BOOST_CHECK(f3.num_parameters() == 3); + BOOST_CHECK(f3[0] == 0); + BOOST_CHECK(f3[1] == 1); + BOOST_CHECK(f3[2] == 2); + + std::vector v{0, 1, 2}; + One_critical_filtration f4(v); + BOOST_CHECK(f4.size() == 3); + BOOST_CHECK(f4.num_parameters() == 3); + BOOST_CHECK(f4[0] == 0); + BOOST_CHECK(f4[1] == 1); + BOOST_CHECK(f4[2] == 2); + + One_critical_filtration f5(v.begin(), v.end()); + BOOST_CHECK(f5.size() == 3); + BOOST_CHECK(f5.num_parameters() == 3); + BOOST_CHECK(f5[0] == 0); + BOOST_CHECK(f5[1] == 1); + BOOST_CHECK(f5[2] == 2); + + f = f5; + BOOST_CHECK(f.size() == 3); + BOOST_CHECK(f.num_parameters() == 3); + BOOST_CHECK(f[0] == 0); + BOOST_CHECK(f[1] == 1); + BOOST_CHECK(f[2] == 2); +} + +BOOST_AUTO_TEST_CASE_TEMPLATE(one_critical_filtration_utilities, T, list_of_tested_variants) +{ + One_critical_filtration f({0, 1, 2}); + bool test = std::is_same_v; + BOOST_CHECK(test); + + One_critical_filtration f2 = f.template as_type(); + test = std::is_same_v; + BOOST_CHECK(test); + BOOST_CHECK(f2.size() == 3); + BOOST_CHECK(f2.num_parameters() == 3); + BOOST_CHECK(f2[0] == 0.); + BOOST_CHECK(f2[1] == 1.); + BOOST_CHECK(f2[2] == 2.); + + BOOST_CHECK(!f.is_plus_inf()); + BOOST_CHECK(!f.is_minus_inf()); + BOOST_CHECK(!f.is_nan()); + BOOST_CHECK(f.is_finite()); + + One_critical_filtration f3; + BOOST_CHECK(!f3.is_plus_inf()); + BOOST_CHECK(f3.is_minus_inf()); + BOOST_CHECK(!f3.is_nan()); + BOOST_CHECK(!f3.is_finite()); + + //{-inf, -inf, -inf} is considered finite as the user is supposed to updates the values to something else + //the idea is just to reserve space and to give the possibility to use `f4[i] =` + //if the value should really be -inf, use `f4(1)` or `f4 = minus_inf()` instead. + One_critical_filtration f4(3); + BOOST_CHECK(!f4.is_plus_inf()); + BOOST_CHECK(!f4.is_minus_inf()); + BOOST_CHECK(!f4.is_nan()); + BOOST_CHECK(f4.is_finite()); + + One_critical_filtration f5(1); + BOOST_CHECK(!f5.is_plus_inf()); + BOOST_CHECK(f5.is_minus_inf()); + BOOST_CHECK(!f5.is_nan()); + BOOST_CHECK(!f5.is_finite()); + + auto f6 = One_critical_filtration::nan(); + BOOST_CHECK(!f6.is_plus_inf()); + BOOST_CHECK(!f6.is_minus_inf()); + BOOST_CHECK(f6.is_nan()); + BOOST_CHECK(!f6.is_finite()); + + auto f7 = One_critical_filtration::inf(); + BOOST_CHECK(f7.is_plus_inf()); + BOOST_CHECK(!f7.is_minus_inf()); + BOOST_CHECK(!f7.is_nan()); + BOOST_CHECK(!f7.is_finite()); +} + +BOOST_AUTO_TEST_CASE_TEMPLATE(one_critical_filtration_comparators, T, list_of_tested_variants) +{ + One_critical_filtration f1({0, 1, 2}); + One_critical_filtration f2({-1, 0, 1}); + One_critical_filtration f3({0, 2, 3}); + One_critical_filtration f4({5, -1, 2}); + + BOOST_CHECK(!(f1 < f1)); + BOOST_CHECK(!(f1 < f2)); + BOOST_CHECK(f1 < f3); + BOOST_CHECK(!(f1 < f4)); + BOOST_CHECK(!(f1 < One_critical_filtration::nan())); + BOOST_CHECK(f1 < One_critical_filtration::inf()); + BOOST_CHECK(!(f1 < One_critical_filtration::minus_inf())); + + BOOST_CHECK(f1 <= f1); + BOOST_CHECK(!(f1 <= f2)); + BOOST_CHECK(f1 <= f3); + BOOST_CHECK(!(f1 <= f4)); + BOOST_CHECK(!(f1 <= One_critical_filtration::nan())); + BOOST_CHECK(f1 <= One_critical_filtration::inf()); + BOOST_CHECK(!(f1 <= One_critical_filtration::minus_inf())); + + BOOST_CHECK(!(f1 > f1)); + BOOST_CHECK(f1 > f2); + BOOST_CHECK(!(f1 > f3)); + BOOST_CHECK(!(f1 > f4)); + BOOST_CHECK(!(f1 > One_critical_filtration::nan())); + BOOST_CHECK(!(f1 > One_critical_filtration::inf())); + BOOST_CHECK(f1 > One_critical_filtration::minus_inf()); + + BOOST_CHECK(f1 >= f1); + BOOST_CHECK(f1 >= f2); + BOOST_CHECK(!(f1 >= f3)); + BOOST_CHECK(!(f1 >= f4)); + BOOST_CHECK(!(f1 >= One_critical_filtration::nan())); + BOOST_CHECK(!(f1 >= One_critical_filtration::inf())); + BOOST_CHECK(f1 >= One_critical_filtration::minus_inf()); + + BOOST_CHECK(f1 == f1); + BOOST_CHECK(!(f1 == f2)); + BOOST_CHECK(!(f1 == f3)); + BOOST_CHECK(!(f1 == f4)); + BOOST_CHECK(!(f1 == One_critical_filtration::nan())); + BOOST_CHECK(!(f1 == One_critical_filtration::inf())); + BOOST_CHECK(!(f1 == One_critical_filtration::minus_inf())); + + BOOST_CHECK(!(f1 != f1)); + BOOST_CHECK(f1 != f2); + BOOST_CHECK(f1 != f3); + BOOST_CHECK(f1 != f4); + BOOST_CHECK(f1 != One_critical_filtration::nan()); + BOOST_CHECK(f1 != One_critical_filtration::inf()); + BOOST_CHECK(f1 != One_critical_filtration::minus_inf()); +} + +BOOST_AUTO_TEST_CASE_TEMPLATE(one_critical_filtration_operators, T, list_of_tested_variants) +{; + One_critical_filtration f({-10, 0, 1}); + One_critical_filtration f2({5, 2, -1}); + One_critical_filtration f3( + {-One_critical_filtration::T_inf, One_critical_filtration::T_inf, -One_critical_filtration::T_inf}); + One_critical_filtration f4( + {One_critical_filtration::T_inf, -One_critical_filtration::T_inf, One_critical_filtration::T_inf}); + + auto res = -f; + BOOST_CHECK_EQUAL(res[0], 10); + BOOST_CHECK_EQUAL(res[1], 0); + BOOST_CHECK_EQUAL(res[2], -1); + BOOST_CHECK((-One_critical_filtration::inf()).is_minus_inf()); + BOOST_CHECK((-One_critical_filtration::minus_inf()).is_plus_inf()); + BOOST_CHECK((-One_critical_filtration::nan()).is_nan()); + + res = f - f2; + BOOST_CHECK_EQUAL(res[0], -15); + BOOST_CHECK_EQUAL(res[1], -2); + BOOST_CHECK_EQUAL(res[2], 2); + + res = f - f3; + BOOST_CHECK_EQUAL(res[0], f4[0]); + BOOST_CHECK_EQUAL(res[1], f4[1]); + BOOST_CHECK_EQUAL(res[2], f4[2]); + + res = f3 - f; + BOOST_CHECK_EQUAL(res[0], f3[0]); + BOOST_CHECK_EQUAL(res[1], f3[1]); + BOOST_CHECK_EQUAL(res[2], f3[2]); + + res = 5 - f; + BOOST_CHECK_EQUAL(res[0], 15); + BOOST_CHECK_EQUAL(res[1], 5); + BOOST_CHECK_EQUAL(res[2], 4); + + res = f - 5; + BOOST_CHECK_EQUAL(res[0], -15); + BOOST_CHECK_EQUAL(res[1], -5); + BOOST_CHECK_EQUAL(res[2], -4); + + BOOST_CHECK((f - One_critical_filtration::inf()).is_minus_inf()); + BOOST_CHECK((One_critical_filtration::inf() - f).is_plus_inf()); + BOOST_CHECK((f - One_critical_filtration::minus_inf()).is_plus_inf()); + BOOST_CHECK((One_critical_filtration::minus_inf() - f).is_minus_inf()); + BOOST_CHECK((f - One_critical_filtration::nan()).is_nan()); + BOOST_CHECK((One_critical_filtration::nan() - f).is_nan()); + + res = f3 - f3; + BOOST_CHECK(res.is_nan()); + res = f3 - f4; + BOOST_CHECK_EQUAL(res[0], f3[0]); + BOOST_CHECK_EQUAL(res[1], f3[1]); + BOOST_CHECK_EQUAL(res[2], f3[2]); + + res = f + f2; + BOOST_CHECK_EQUAL(res[0], -5); + BOOST_CHECK_EQUAL(res[1], 2); + BOOST_CHECK_EQUAL(res[2], 0); + + res = f + f3; + BOOST_CHECK_EQUAL(res[0], f3[0]); + BOOST_CHECK_EQUAL(res[1], f3[1]); + BOOST_CHECK_EQUAL(res[2], f3[2]); + + res = 5 + f; + BOOST_CHECK_EQUAL(res[0], -5); + BOOST_CHECK_EQUAL(res[1], 5); + BOOST_CHECK_EQUAL(res[2], 6); + + res = f + 5; + BOOST_CHECK_EQUAL(res[0], -5); + BOOST_CHECK_EQUAL(res[1], 5); + BOOST_CHECK_EQUAL(res[2], 6); + + BOOST_CHECK((f + One_critical_filtration::inf()).is_plus_inf()); + BOOST_CHECK((One_critical_filtration::inf() + f).is_plus_inf()); + BOOST_CHECK((f + One_critical_filtration::minus_inf()).is_minus_inf()); + BOOST_CHECK((One_critical_filtration::minus_inf() + f).is_minus_inf()); + BOOST_CHECK((f + One_critical_filtration::nan()).is_nan()); + BOOST_CHECK((One_critical_filtration::nan() + f).is_nan()); + + res = f3 + f4; + BOOST_CHECK(res.is_nan()); + res = f3 + f3; + BOOST_CHECK_EQUAL(res[0], f3[0]); + BOOST_CHECK_EQUAL(res[1], f3[1]); + BOOST_CHECK_EQUAL(res[2], f3[2]); + + res = f * f2; + BOOST_CHECK_EQUAL(res[0], -50); + BOOST_CHECK_EQUAL(res[1], 0); + BOOST_CHECK_EQUAL(res[2], -1); + + res = f * f3; + if constexpr (std::numeric_limits::has_quiet_NaN){ + BOOST_CHECK_EQUAL(res[0], f4[0]); + BOOST_CHECK(std::isnan(res[1])); + BOOST_CHECK_EQUAL(res[2], f3[2]); + } else { + BOOST_CHECK(res.is_nan()); + } + + res = 5 * f; + BOOST_CHECK_EQUAL(res[0], -50); + BOOST_CHECK_EQUAL(res[1], 0); + BOOST_CHECK_EQUAL(res[2], 5); + + res = f * 5; + BOOST_CHECK_EQUAL(res[0], -50); + BOOST_CHECK_EQUAL(res[1], 0); + BOOST_CHECK_EQUAL(res[2], 5); + + res = f * One_critical_filtration::inf(); + if constexpr (std::numeric_limits::has_quiet_NaN){ + BOOST_CHECK_EQUAL(res[0], -One_critical_filtration::T_inf); + BOOST_CHECK(std::isnan(res[1])); + BOOST_CHECK_EQUAL(res[2], One_critical_filtration::T_inf); + } else { + BOOST_CHECK(res.is_nan()); + } + res = One_critical_filtration::inf() * f; + if constexpr (std::numeric_limits::has_quiet_NaN){ + BOOST_CHECK_EQUAL(res[0], -One_critical_filtration::T_inf); + BOOST_CHECK(std::isnan(res[1])); + BOOST_CHECK_EQUAL(res[2], One_critical_filtration::T_inf); + } else { + BOOST_CHECK(res.is_nan()); + } + res = f * One_critical_filtration::minus_inf(); + if constexpr (std::numeric_limits::has_quiet_NaN){ + BOOST_CHECK_EQUAL(res[0], One_critical_filtration::T_inf); + BOOST_CHECK(std::isnan(res[1])); + BOOST_CHECK_EQUAL(res[2], -One_critical_filtration::T_inf); + } else { + BOOST_CHECK(res.is_nan()); + } + res = One_critical_filtration::minus_inf() * f; + if constexpr (std::numeric_limits::has_quiet_NaN){ + BOOST_CHECK_EQUAL(res[0], One_critical_filtration::T_inf); + BOOST_CHECK(std::isnan(res[1])); + BOOST_CHECK_EQUAL(res[2], -One_critical_filtration::T_inf); + } else { + BOOST_CHECK(res.is_nan()); + } + res = f * One_critical_filtration::nan(); + BOOST_CHECK(res.is_nan()); + res = One_critical_filtration::nan() * f; + BOOST_CHECK(res.is_nan()); + + res = f3 * f3; + BOOST_CHECK(res.is_plus_inf()); + res = f3 * f4; + BOOST_CHECK(res.is_minus_inf()); + + res = f / f2; + BOOST_CHECK_EQUAL(res[0], -2); + BOOST_CHECK_EQUAL(res[1], 0); + BOOST_CHECK_EQUAL(res[2], -1); + + res = f / f3; + BOOST_CHECK_EQUAL(res[0], 0); + BOOST_CHECK_EQUAL(res[1], 0); + BOOST_CHECK_EQUAL(res[2], 0); + + res = f3 / f; + if constexpr (std::numeric_limits::has_quiet_NaN){ + BOOST_CHECK_EQUAL(res[0], f4[0]); + BOOST_CHECK(std::isnan(res[1])); + BOOST_CHECK_EQUAL(res[2], f3[2]); + } else { + BOOST_CHECK(res.is_nan()); + } + + res = 5 / f; + if constexpr (std::numeric_limits::has_quiet_NaN){ + BOOST_CHECK_EQUAL(res[0], -0.5); + BOOST_CHECK(std::isnan(res[1])); + BOOST_CHECK_EQUAL(res[2], 5); + } else { + BOOST_CHECK(res.is_nan()); + } + + res = f / 5; + BOOST_CHECK_EQUAL(res[0], -2); + BOOST_CHECK_EQUAL(res[1], 0); + BOOST_CHECK_EQUAL(res[2], static_cast(1) / static_cast(5)); //to avoid precision error + + res = f / One_critical_filtration::inf(); + BOOST_CHECK_EQUAL(res[0], 0); + BOOST_CHECK_EQUAL(res[1], 0); + BOOST_CHECK_EQUAL(res[2], 0); + res = One_critical_filtration::inf() / f; + if constexpr (std::numeric_limits::has_quiet_NaN){ + BOOST_CHECK_EQUAL(res[0], -One_critical_filtration::T_inf); + BOOST_CHECK(std::isnan(res[1])); + BOOST_CHECK_EQUAL(res[2], One_critical_filtration::T_inf); + } else { + BOOST_CHECK(res.is_nan()); + } + res = f / One_critical_filtration::minus_inf(); + BOOST_CHECK_EQUAL(res[0], 0); + BOOST_CHECK_EQUAL(res[1], 0); + BOOST_CHECK_EQUAL(res[2], 0); + res = One_critical_filtration::minus_inf() / f; + if constexpr (std::numeric_limits::has_quiet_NaN){ + BOOST_CHECK_EQUAL(res[0], One_critical_filtration::T_inf); + BOOST_CHECK(std::isnan(res[1])); + BOOST_CHECK_EQUAL(res[2], -One_critical_filtration::T_inf); + } else { + BOOST_CHECK(res.is_nan()); + } + res = f / One_critical_filtration::nan(); + BOOST_CHECK(res.is_nan()); + res = One_critical_filtration::nan() / f; + BOOST_CHECK(res.is_nan()); + + res = f3 / f3; + BOOST_CHECK(res.is_nan()); + res = f3 / f4; + BOOST_CHECK(res.is_nan()); + res = f / One_critical_filtration({0, 0, 0}); + BOOST_CHECK(res.is_nan()); +} + +BOOST_AUTO_TEST_CASE_TEMPLATE(one_critical_filtration_modifiers, T, list_of_tested_variants) +{ + One_critical_filtration f({0, 1, 2}); + BOOST_CHECK_EQUAL(f[0], 0); + BOOST_CHECK_EQUAL(f[1], 1); + BOOST_CHECK_EQUAL(f[2], 2); + + f.push_to_least_common_upper_bound({-1, 5, 6}); + BOOST_CHECK_EQUAL(f[0], 0); + BOOST_CHECK_EQUAL(f[1], 5); + BOOST_CHECK_EQUAL(f[2], 6); + + f.push_to_least_common_upper_bound({-1, -5, -6}); + BOOST_CHECK_EQUAL(f[0], 0); + BOOST_CHECK_EQUAL(f[1], 5); + BOOST_CHECK_EQUAL(f[2], 6); + + f.push_to_least_common_upper_bound(One_critical_filtration::minus_inf()); + BOOST_CHECK_EQUAL(f[0], 0); + BOOST_CHECK_EQUAL(f[1], 5); + BOOST_CHECK_EQUAL(f[2], 6); + + f.push_to_least_common_upper_bound(One_critical_filtration::inf()); + BOOST_CHECK(f.is_plus_inf()); + + f.push_to_least_common_upper_bound(One_critical_filtration::nan()); + BOOST_CHECK(f.is_plus_inf()); + + f.pull_to_greatest_common_lower_bound({-1, 5, 6}); + BOOST_CHECK_EQUAL(f[0], -1); + BOOST_CHECK_EQUAL(f[1], 5); + BOOST_CHECK_EQUAL(f[2], 6); + + f.pull_to_greatest_common_lower_bound({1, 8, 9}); + BOOST_CHECK_EQUAL(f[0], -1); + BOOST_CHECK_EQUAL(f[1], 5); + BOOST_CHECK_EQUAL(f[2], 6); + + f.pull_to_greatest_common_lower_bound(One_critical_filtration::inf()); + BOOST_CHECK_EQUAL(f[0], -1); + BOOST_CHECK_EQUAL(f[1], 5); + BOOST_CHECK_EQUAL(f[2], 6); + + f.pull_to_greatest_common_lower_bound(One_critical_filtration::minus_inf()); + BOOST_CHECK(f.is_minus_inf()); + + f.pull_to_greatest_common_lower_bound(One_critical_filtration::nan()); + BOOST_CHECK(f.is_minus_inf()); + + std::vector > grid = {{0, 2, 4, 8}, {0, 3, 6, 9}, {0, 4, 8, 16}}; + + f.push_to_least_common_upper_bound({1, 7, 5}); + f.project_onto_grid(grid, true); + BOOST_CHECK_EQUAL(f[0], 1); + BOOST_CHECK_EQUAL(f[1], 3); + BOOST_CHECK_EQUAL(f[2], 2); + + f.push_to_least_common_upper_bound({1, 7, 5}); + f.project_onto_grid(grid, false); + BOOST_CHECK_EQUAL(f[0], 2); + BOOST_CHECK_EQUAL(f[1], 9); + BOOST_CHECK_EQUAL(f[2], 8); +} + +BOOST_AUTO_TEST_CASE_TEMPLATE(one_critical_filtration_friends, T, list_of_tested_variants) +{ + One_critical_filtration f({0, 1, 2}); + + BOOST_CHECK_EQUAL(compute_linear_projection(f, {2,3,5,9}), 13); + BOOST_CHECK_EQUAL(compute_norm(f), static_cast(std::sqrt(T(5)))); + BOOST_CHECK_EQUAL(compute_euclidean_distance_to(f, {2,3,5}), static_cast(std::sqrt(T(17)))); + + f = {1, 7, 5}; + + std::vector > grid = {{0, 2, 4, 8}, {0, 3, 6, 9}, {0, 4, 8, 16}}; + auto res = compute_coordinates_in_grid(f, grid); + BOOST_CHECK_EQUAL(res[0], 1); + BOOST_CHECK_EQUAL(res[1], 3); + BOOST_CHECK_EQUAL(res[2], 2); + + res = evaluate_coordinates_in_grid(res, grid); + BOOST_CHECK_EQUAL(res[0], 2); + BOOST_CHECK_EQUAL(res[1], 9); + BOOST_CHECK_EQUAL(res[2], 8); +} + +BOOST_AUTO_TEST_CASE_TEMPLATE(one_critical_filtration_numerical_limits, T, list_of_tested_variants) +{ + BOOST_CHECK(std::numeric_limits >::has_infinity); + BOOST_CHECK(std::numeric_limits >::infinity().is_plus_inf()); + BOOST_CHECK(std::numeric_limits >::quiet_NaN().is_nan()); + BOOST_CHECK_THROW(std::numeric_limits >::max(), std::logic_error); + auto max = std::numeric_limits >::max(3); + BOOST_CHECK_EQUAL(max[0], std::numeric_limits::max()); + BOOST_CHECK_EQUAL(max[1], std::numeric_limits::max()); + BOOST_CHECK_EQUAL(max[2], std::numeric_limits::max()); +} + diff --git a/src/Multi_persistence/doc/COPYRIGHT b/src/Multi_persistence/doc/COPYRIGHT new file mode 100644 index 0000000000..70cac8b762 --- /dev/null +++ b/src/Multi_persistence/doc/COPYRIGHT @@ -0,0 +1,12 @@ +The files of this directory are 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): Hannah Schreiber, David Loiseaux + +Copyright (C) 2024 Inria + +This gives everyone the freedoms to use openFrameworks in any context: +commercial or non-commercial, public or private, open or closed source. + +You should have received a copy of the MIT License along with this program. +If not, see https://opensource.org/licenses/MIT. diff --git a/src/Multi_persistence/doc/Intro_multi_persistence.h b/src/Multi_persistence/doc/Intro_multi_persistence.h new file mode 100644 index 0000000000..a23bed4382 --- /dev/null +++ b/src/Multi_persistence/doc/Intro_multi_persistence.h @@ -0,0 +1,31 @@ +/* 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): Hannah Schreiber + * + * Copyright (C) 2024 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#ifndef DOC_MULTI_PERSISTENCE_INTRO_H_ +#define DOC_MULTI_PERSISTENCE_INTRO_H_ + +// needs namespace for Doxygen to link on classes +namespace Gudhi { +namespace multi_persistence { + +/** \defgroup multi_persistence Multi-parameter Persistence + * @{ + * \author David Loiseaux + * + * \section multiparameterintro Multi-parameter Persistence + * + * To come: a Multi-parameter Persistence module, with, for the beginnings, an interface mostly in Python. + * + * @} + */ +} // namespace persistence_fields +} // namespace Gudhi + +#endif // DOC_MULTI_PERSISTENCE_INTRO_H_ diff --git a/src/Multi_persistence/doc/multidiagrams.png b/src/Multi_persistence/doc/multidiagrams.png new file mode 100644 index 0000000000..2cf004e46e Binary files /dev/null and b/src/Multi_persistence/doc/multidiagrams.png differ diff --git a/src/Multi_persistence/include/gudhi/Multi_persistence/Box.h b/src/Multi_persistence/include/gudhi/Multi_persistence/Box.h new file mode 100644 index 0000000000..2f90f4d99a --- /dev/null +++ b/src/Multi_persistence/include/gudhi/Multi_persistence/Box.h @@ -0,0 +1,169 @@ +/* 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): David Loiseaux + * + * Copyright (C) 2023 Inria + * + * Modification(s): + * - 2024/08 Hannah Schreiber: doc + * - YYYY/MM Author: Description of the modification + */ + +/** + * @file Box.h + * @author David Loiseaux + * @brief Contains the @ref Gudhi::multi_persistence::Box class. + */ + +#ifndef BOX_H_INCLUDED +#define BOX_H_INCLUDED + +#include //std::ostream + +#include +#include + +namespace Gudhi::multi_persistence { + +/** + * @class Box Box.h gudhi/Multi_persistence/Box.h + * @ingroup multi_persistence + * + * @brief Simple box in \f$\mathbb R^n\f$ defined by two diametrically opposite corners. + * + * @tparam T Type of the coordinates of the Box. Has to follow the conditions of the template parameter of + * @ref One_critical_filtration "". + */ +template +class Box { + public: + using Point = Gudhi::multi_filtration::One_critical_filtration; /**< Type of a point in \f$\mathbb R^n\f$. */ + + /** + * @brief Default constructor. Constructs a trivial box with corners at minus infinity. + */ + Box() {} + + /** + * @brief Constructs a box from the two given corners. Assumes that \f$ lowerCorner \le @p upperCorner \f$ and + * if both are finite values, they have the same dimension. + * + * @param lowerCorner First corner of the box. Has to be smaller than `upperCorner`. + * @param upperCorner Second corner of the box. Has to be greater than `lowerCorner`. + */ + Box(const Point &lowerCorner, const Point &upperCorner) : lowerCorner_(lowerCorner), upperCorner_(upperCorner) { + GUDHI_CHECK(lowerCorner.size() == upperCorner.size() && lowerCorner <= upperCorner, "This box is trivial !"); + } + + /** + * @brief Constructs a box from the two given corners. Assumes that \f$ box.first \le @p box.second \f$ and + * if both are finite values, they have the same dimension. + * + * @param box Pair of corners defining the wished box. + */ + Box(const std::pair &box) : Box(box.first, box.second) {} + + /** + * @brief Returns the lowest of both defining corners. + */ + const Point &get_lower_corner() const { return lowerCorner_; } + + /** + * @brief Returns the lowest of both defining corners. + */ + Point &get_lower_corner() { return lowerCorner_; } + + /** + * @brief Returns the greatest of both defining corners. + */ + Point &get_upper_corner() { return upperCorner_; } + + /** + * @brief Returns the greatest of both defining corners. + */ + const Point &get_upper_corner() const { return upperCorner_; } + + /** + * @brief Returns a pair of const references to both defining corners. + */ + std::pair get_bounding_corners() const { return {lowerCorner_, upperCorner_}; } + + /** + * @brief Returns a pair of references to both defining corners. + */ + std::pair get_bounding_corners() { return {lowerCorner_, upperCorner_}; } + + /** + * @brief Returns true if and only if one of the following is true: + * - one of the corners is empty + * - one of the corners has value NaN + * - both corners have value infinity + * - both corners have value minus infinity + * - both corners are finite but don't have the same dimension. + */ + bool is_trivial() const { + return lowerCorner_.empty() || upperCorner_.empty() || lowerCorner_.is_nan() || upperCorner_.is_nan() || + (lowerCorner_.is_plus_inf() && upperCorner_.is_plus_inf()) || + (lowerCorner_.is_minus_inf() && upperCorner_.is_minus_inf()) || + (lowerCorner_.is_finite() && upperCorner_.is_finite() && + lowerCorner_.num_parameters() != upperCorner_.num_parameters()); + } + + /** + * @brief Returns true if and only if the given point is inside the box. + * If the box is not {-infinity, infinity} and the given point is finite, but has not the same dimension + * than the box, the point is considered outside. + */ + bool contains(const Point &point) const { + if (point.is_nan() || is_trivial()) return false; + if (point.is_plus_inf()) return upperCorner_.is_plus_inf(); + if (point.is_minus_inf()) return lowerCorner_.is_minus_inf(); + + if ((lowerCorner_.is_finite() && point.size() != lowerCorner_.size()) || + (upperCorner_.is_finite() && point.size() != upperCorner_.size())) { + // TODO: make it a warning, with future GUDHI_CHECK version? + // std::cerr << "Box and point are not of the same dimension." << std::endl; + return false; + } + + return lowerCorner_ <= point && point <= upperCorner_; + } + + /** + * @brief Returns the dimension of the box. If the box is trivial or both corners are infinite, the dimension is 0. + */ + std::size_t dimension() const { + if (is_trivial()) return 0; + if (lowerCorner_.is_minus_inf() && upperCorner_.is_plus_inf()) return 0; // not so sure what we want to do here + return lowerCorner_.is_finite() ? lowerCorner_.size() : upperCorner_.size(); + } + + /** + * @brief Inflates the box by delta. + * + * @param delta Inflation coefficient. + */ + void inflate(T delta) { + lowerCorner_ -= delta; + upperCorner_ += delta; + } + + /** + * @brief Outstream operator. + */ + friend std::ostream &operator<<(std::ostream &os, const Box &box) { + os << "Box -- Bottom corner : "; + os << box.get_lower_corner(); + os << ", Top corner : "; + os << box.get_upper_corner(); + return os; + } + + private: + Point lowerCorner_; /**< Lowest of defining corners. */ + Point upperCorner_; /**< Greatest of defining corners. */ +}; + +} // namespace Gudhi::multi_persistence + +#endif // BOX_H_INCLUDED diff --git a/src/Multi_persistence/include/gudhi/Multi_persistence/Line.h b/src/Multi_persistence/include/gudhi/Multi_persistence/Line.h new file mode 100644 index 0000000000..378e968366 --- /dev/null +++ b/src/Multi_persistence/include/gudhi/Multi_persistence/Line.h @@ -0,0 +1,281 @@ +/* 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): David Loiseaux + * + * Copyright (C) 2023 Inria + * + * Modification(s): + * - 2024/08 Hannah Schreiber: doc + * - YYYY/MM Author: Description of the modification + */ + +/** + * @file Line.h + * @author David Loiseaux + * @brief Contains the @ref Gudhi::multi_persistence::Line class. + */ + +#ifndef LINE_FILTRATION_TRANSLATION_H_INCLUDED +#define LINE_FILTRATION_TRANSLATION_H_INCLUDED + +#include +#include +#include +#include + +#include +#include +#include +#include + +namespace Gudhi::multi_persistence { + +/** + * @class Line Line.h gudhi/Multi_persistence/Line.h + * @ingroup multi_persistence + * + * @brief A line in \f$\mathbb R^n\f$, with some helpers to project points on it. + * + * @tparam T Type of the coordinate values. Has to follow the conditions of the template parameter of + * @ref One_critical_filtration "". + */ +template +class Line { + public: + /** + * @brief Coordinates in \f$\mathbb R^n\f$. + */ + using Point = Gudhi::multi_filtration::One_critical_filtration; + /** + * @brief Set of coordinates in \f$\mathbb R^n\f$. + */ + using K_critical_point = Gudhi::multi_filtration::Multi_critical_filtration; + + /** + * @brief Default constructor. Sets the number of coordinates to 0. + */ + Line() : basePoint_(0), direction_(0) {} // has to be explicitly set to 0, otherwise becomes -inf + /** + * @brief Constructs a line going through the given point with slope 1. + * + * @param x A point of the line. + */ + Line(const Point &x) : basePoint_(x), direction_(0) {} // default direction + /** + * @brief Constructs a line going through the given point with slope 1. + * + * @param x A point of the line. Will be moved. + */ + Line(Point &&x) : basePoint_(std::move(x)), direction_(0) {} // default direction + /** + * @brief Constructs a line going through the given point in the direction of the given vector. + * If the vector has no coordinates, the slope is assumed to be 1. + * Otherwise, the vector has to be non trivial and all its coordinates have to be positive. + * + * @param x A point of the line. + * @param vector Direction of the line. Positive and non trivial. + */ + Line(const Point &x, const Point &vector) : basePoint_(x), direction_(vector) { check_direction_(); } + + /** + * @brief Returns the coordinates of the point on the line with "time" parameter `t`. That is, the point \f$ x \f$ + * such that \f$ x[i] = base\_point[i] + t \times direction[i] \f$ for all \f$ i \in [0, n - 1] \f$ with \f$ n \f$ + * the number of coordinates. + */ + Point operator[](T t) const { + GUDHI_CHECK(direction_.empty() || direction_.size() == basePoint_.size(), + "Direction and base point do not have the same dimension."); + + if constexpr (std::numeric_limits::has_quiet_NaN){ //to avoid windows error + if (std::isnan(t)) return Point::nan(); + } + if (t == Point::T_inf) return Point::inf(); + if (t == -Point::T_inf) return Point::minus_inf(); + + Point x(basePoint_.size()); + + if (direction_.size() > 0) { + for (std::size_t i = 0; i < x.size(); i++) x[i] = basePoint_[i] + t * direction_[i]; + } else + for (std::size_t i = 0; i < x.size(); i++) x[i] = basePoint_[i] + t; + + return x; + } + + /** + * @brief Translates the given line in the given direction. + */ + friend Line &operator+=(Line &to_translate, const Point &v) { + to_translate.basePoint_ += v; + return to_translate; + } + + /** + * @brief Returns a reference to the current base point of the line. + */ + Point &base_point() { return basePoint_; } + /** + * @brief Returns a const reference to the current base point of the line. + */ + const Point &base_point() const { return basePoint_; } + + /** + * @brief Returns a reference to the direction vector of the line. + */ + Point &direction() { return direction_; } + /** + * @brief Returns a const reference to the direction vector of the line. + */ + const Point &direction() const { return direction_; } + + // TODO: factorize forward and backward version by adding a `co` to One_critical_filtration? + // Could make problems with One_critical_filtration being the type of basePoint_ and direction_ + + /** + * @brief Computes the "time" parameter \f$ t \f$ of the starting point \f$ p = base\_point + t \times direction \f$ + * of the intersection between the line and the closed positive cone originating at `x`. + * + * @tparam U Type of the time parameter. + * @param x Origin of the closed positive cone. + */ + template + U compute_forward_intersection(const Point &x) const { + GUDHI_CHECK(direction_.empty() || direction_.size() == x.size(), "x has not as many parameters as the line."); + + constexpr const U inf = + std::numeric_limits::has_infinity ? std::numeric_limits::infinity() : std::numeric_limits::max(); + if (x.is_plus_inf() || x.is_nan()) return inf; + if (x.is_minus_inf()) return -inf; + U t = -inf; + if (direction_.size()) { + for (std::size_t i = 0; i < x.size(); i++) { + if (direction_[i] == 0) { + if (x[i] > basePoint_[i]) return inf; + } else { + t = std::max(t, (static_cast(x[i]) - static_cast(basePoint_[i])) / static_cast((direction_[i]))); + } + } + } else { + for (std::size_t i = 0; i < x.size(); i++) t = std::max(t, static_cast(x[i]) - static_cast(basePoint_[i])); + } + + return t; + } + + /** + * @brief Computes the "time" parameter \f$ t \f$ of the starting point \f$ p = base\_point + t \times direction \f$ + * of the intersection between the line and the union of closed positive cones originating at the points in `x`. + * + * @tparam U Type of the time parameter. + * @param x Set of origins for the closed positive cones. + */ + template + U compute_forward_intersection(const K_critical_point &x) const { + constexpr const U inf = + std::numeric_limits::has_infinity ? std::numeric_limits::infinity() : std::numeric_limits::max(); + if (x.is_plus_inf() || x.is_nan()) return inf; + if (x.is_minus_inf()) return -inf; + U t = inf; + for (const auto &y : x) { + t = std::min(t, compute_forward_intersection(y)); + } + return t; + } + + /** + * @brief Computes the "time" parameter \f$ t \f$ of the starting point \f$ p = base\_point + t \times direction \f$ + * of the intersection between the line and the open negative cone originating at `x`. + * + * @tparam U Type of the time parameter. + * @param x Origin of the open negative cone. + */ + template + U compute_backward_intersection(const Point &x) const { + constexpr const U inf = + std::numeric_limits::has_infinity ? std::numeric_limits::infinity() : std::numeric_limits::max(); + if (x.is_plus_inf()) return inf; + if (x.is_minus_inf() || x.is_nan()) return -inf; + U t = inf; + + if (direction_.size()) { + for (std::size_t i = 0; i < x.size(); i++) { + if (direction_[i] == 0) { + if (x[i] <= basePoint_[i]) return -inf; + } else { + t = std::min(t, (static_cast(x[i]) - static_cast(basePoint_[i])) / static_cast(direction_[i])); + } + } + } else { + for (std::size_t i = 0; i < x.size(); i++) t = std::min(t, static_cast(x[i] - basePoint_[i])); + } + return t; + } + + /** + * @brief Computes the "time" parameter \f$ t \f$ of the starting point \f$ p = base\_point + t \times direction \f$ + * of the intersection between the line and the union of open negative cones originating at the points in `x`. + * + * @tparam U Type of the time parameter. + * @param x Set of origins for the open negative cones. + */ + template + U compute_backward_intersection(const K_critical_point &x) const { + constexpr const U inf = + std::numeric_limits::has_infinity ? std::numeric_limits::infinity() : std::numeric_limits::max(); + if (x.is_plus_inf()) return inf; + if (x.is_minus_inf() || x.is_nan()) return -inf; + U t = -inf; + for (const auto &y : x) { + t = std::max(t, compute_backward_intersection(y)); + } + return t; + } + + /** + * @brief Given a box, returns "time" parameter of the intersection of this box and the line. + * + * @param box Box to intersect. + * @return A pair representing the two bounding points of the intersection, such that the first element is the + * smallest of the two. If the box and the line do not intersect or the box is trivial, returns the pair {inf, -inf}. + */ + std::pair get_bounds(const Box &box) const { + if (box.is_trivial()) return {Point::T_inf, -Point::T_inf}; + + T bottom = compute_forward_intersection(box.get_lower_corner()); + T top = compute_backward_intersection(box.get_upper_corner()); + + if (bottom > top) return {Point::T_inf, -Point::T_inf}; // no intersection + + return {bottom, top}; + } + + private: + Point basePoint_; /**< Any point on the line. */ + Point direction_; /**< Direction of the line. */ + + /** + * @brief Checks that the arguments define a correct and positively slopped line. + */ + void check_direction_() const { + if (direction_.size() == 0) return; // default slope + + bool is_trivial = true; + for (T v : direction_) { + if (v) { + is_trivial = false; + } + if (v < 0) { + throw std::invalid_argument("Direction should have positive entries."); + } + } + if (is_trivial) { + throw std::invalid_argument("Direction should have at least one non-trivial entry."); + } + if (direction_.size() != basePoint_.size()) + throw std::invalid_argument("The dimensions of base point and direction are not equal."); + } +}; + +} // namespace Gudhi::multi_persistence + +#endif // LINE_FILTRATION_TRANSLATION_H_INCLUDED diff --git a/src/Multi_persistence/include/gudhi/Simplex_tree_multi.h b/src/Multi_persistence/include/gudhi/Simplex_tree_multi.h new file mode 100644 index 0000000000..fd8a27fa53 --- /dev/null +++ b/src/Multi_persistence/include/gudhi/Simplex_tree_multi.h @@ -0,0 +1,131 @@ +/* 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): David + * Loiseaux, Hannah Schreiber + * + * Copyright (C) 2023 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ +#ifndef SIMPLEX_TREE_MULTI_H_ +#define SIMPLEX_TREE_MULTI_H_ + +#include +#include + +namespace Gudhi::multi_persistence { + +/** Model of SimplexTreeOptions, with a multiparameter filtration. + * \ingroup multi_persistence + * */ +template +struct Simplex_tree_options_multidimensional_filtration { +public: + typedef linear_indexing_tag Indexing_tag; + typedef int Vertex_handle; + using Filtration_value = Filtration; + typedef typename Filtration::value_type value_type; + typedef std::uint32_t Simplex_key; + static const bool store_key = true; + static const bool store_filtration = true; + static const bool contiguous_vertices = false; + static const bool link_nodes_by_label = false; + static const bool stable_simplex_handles = false; + static const bool is_multi_parameter = true; +}; + +/** + * \brief Turns a 1-parameter simplextree into a multiparameter simplextree, + * and keeps the 1-filtration in the 1st axis. + * Default values can be specified. + * \ingroup multi_persistence + * \tparam simplextree_std A non-multi simplextree + * \tparam simplextree_multi A multi simplextree + * \param st Simplextree to copy + * \param st_multi Multiparameter simplextree container to fill. + * \param default_values If given, this vector is assume to be of size `num_parameters-1` and contains the default + * values of axes `1` to `num_parameters`. + * */ +template +void multify(simplextree_std &st, simplextree_multi &st_multi, + const int num_parameters, + const typename simplextree_multi::Options::Filtration_value + &default_values = {}) { + typename simplextree_multi::Options::Filtration_value f(num_parameters); + static_assert( + !simplextree_std::Options::is_multi_parameter && + simplextree_multi::Options::is_multi_parameter, + "Can only convert non-multiparameter to multiparameter simplextree."); + unsigned int num_default_values; + if constexpr (simplextree_multi::Options::Filtration_value:: + is_multi_critical) { + num_default_values = default_values[0].size(); + } else { + num_default_values = default_values.size(); + } + for (auto i = 0u; i < std::min(num_default_values, + static_cast(num_parameters - 1)); + i++) + if constexpr (simplextree_multi::Options::Filtration_value:: + is_multi_critical) { + f[0][i + 1] = default_values[0][i]; + } else { + f[i + 1] = default_values[i]; + } + + std::vector simplex; + simplex.reserve(st.dimension() + 1); + for (auto &simplex_handle : st.complex_simplex_range()) { + simplex.clear(); + for (auto vertex : st.simplex_vertex_range(simplex_handle)) + simplex.push_back(vertex); + + if (num_parameters > 0) { + if constexpr (simplextree_multi::Options::Filtration_value:: + is_multi_critical) { + f[0][0] = st.filtration(simplex_handle); + } else { + f[0] = st.filtration(simplex_handle); + } + } + st_multi.insert_simplex(simplex, f); + } + st_multi.set_number_of_parameters(num_parameters); +} + +/** + * \brief Turns a multiparameter-parameter simplextree into a 1-parameter + * simplextree. + * \ingroup multi_persistence + * \tparam simplextree_std A non-multi simplextree + * \tparam simplextree_multi A multi simplextree + * \param st Simplextree to fill. + * \param st_multi Multiparameter simplextree to convert into a 1 parameter simplex tree. + * \param dimension The filtration parameter to put into the 1 parameter simplextree. + * */ +template +void flatten(simplextree_std &st, simplextree_multi &st_multi, + const int dimension = 0) { + static_assert( + !simplextree_std::Options::is_multi_parameter && + simplextree_multi::Options::is_multi_parameter, + "Can only convert multiparameter to non-multiparameter simplextree."); + for (const auto &simplex_handle : st_multi.complex_simplex_range()) { + std::vector simplex; + typename simplextree_multi::Options::value_type f; + for (auto vertex : st_multi.simplex_vertex_range(simplex_handle)) + simplex.push_back(vertex); + if constexpr (simplextree_multi::Filtration_value::is_multi_critical) { + f = dimension >= 0 ? st_multi.filtration(simplex_handle)[0][dimension] + : 0; + } else { + f = dimension >= 0 ? st_multi.filtration(simplex_handle)[dimension] : 0; + } + st.insert_simplex(simplex, f); + } +} + +} // namespace Gudhi::multi_persistence + +#endif // SIMPLEX_TREE_MULTI_H_ diff --git a/src/Multi_persistence/test/CMakeLists.txt b/src/Multi_persistence/test/CMakeLists.txt new file mode 100644 index 0000000000..db3e67c346 --- /dev/null +++ b/src/Multi_persistence/test/CMakeLists.txt @@ -0,0 +1,7 @@ +include(GUDHI_boost_test) + +add_executable_with_targets(Multi_persistence_box_unit_test multipersistence_box_unit_test.cpp TBB::tbb) +gudhi_add_boost_test(Multi_persistence_box_unit_test) + +add_executable_with_targets(Multi_persistence_line_unit_test multipersistence_line_unit_test.cpp TBB::tbb) +gudhi_add_boost_test(Multi_persistence_line_unit_test) diff --git a/src/Multi_persistence/test/multipersistence_box_unit_test.cpp b/src/Multi_persistence/test/multipersistence_box_unit_test.cpp new file mode 100644 index 0000000000..353df8aebe --- /dev/null +++ b/src/Multi_persistence/test/multipersistence_box_unit_test.cpp @@ -0,0 +1,140 @@ +/* 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): Hannah Schreiber + * + * Copyright (C) 2024 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#define BOOST_TEST_DYN_LINK +#define BOOST_TEST_MODULE "multi_persistence" +#include +#include +#include + +#include + +using Gudhi::multi_persistence::Box; + +typedef boost::mpl::list list_of_tested_variants; + +BOOST_AUTO_TEST_CASE_TEMPLATE(box_constructors, T, list_of_tested_variants) +{ + Box b; + BOOST_CHECK(b.is_trivial()); + auto bottom = b.get_lower_corner(); + auto top = b.get_upper_corner(); + BOOST_CHECK_EQUAL(bottom.size(), 1); + BOOST_CHECK_EQUAL(top.size(), 1); + BOOST_CHECK(bottom.is_minus_inf()); + BOOST_CHECK(top.is_minus_inf()); + BOOST_CHECK(b.get_bounding_corners().first == bottom); + BOOST_CHECK(b.get_bounding_corners().second == top); + + Box b2({0,1,2}, {4,5,6}); + BOOST_CHECK(!b2.is_trivial()); + bottom = b2.get_lower_corner(); + top = b2.get_upper_corner(); + BOOST_CHECK_EQUAL(bottom.size(), 3); + BOOST_CHECK_EQUAL(top.size(), 3); + BOOST_CHECK_EQUAL(bottom[0], 0); + BOOST_CHECK_EQUAL(bottom[1], 1); + BOOST_CHECK_EQUAL(bottom[2], 2); + BOOST_CHECK_EQUAL(top[0], 4); + BOOST_CHECK_EQUAL(top[1], 5); + BOOST_CHECK_EQUAL(top[2], 6); + BOOST_CHECK(b2.get_bounding_corners().first == bottom); + BOOST_CHECK(b2.get_bounding_corners().second == top); + + Box b3(std::make_pair::Point, typename Box::Point>({0,1,2}, {4,5,6})); + BOOST_CHECK(!b3.is_trivial()); + bottom = b3.get_lower_corner(); + top = b3.get_upper_corner(); + BOOST_CHECK_EQUAL(bottom.size(), 3); + BOOST_CHECK_EQUAL(top.size(), 3); + BOOST_CHECK_EQUAL(bottom[0], 0); + BOOST_CHECK_EQUAL(bottom[1], 1); + BOOST_CHECK_EQUAL(bottom[2], 2); + BOOST_CHECK_EQUAL(top[0], 4); + BOOST_CHECK_EQUAL(top[1], 5); + BOOST_CHECK_EQUAL(top[2], 6); + BOOST_CHECK(b3.get_bounding_corners().first == bottom); + BOOST_CHECK(b3.get_bounding_corners().second == top); +} + +BOOST_AUTO_TEST_CASE_TEMPLATE(box_other, T, list_of_tested_variants) +{ + using P = typename Box::Point; + + Box b({0,1,2}, {4,5,6}); + + BOOST_CHECK(b.contains({1, 2, 3})); + BOOST_CHECK(!b.contains({1, 8, 3})); + BOOST_CHECK(!b.contains(P::nan())); + BOOST_CHECK(!b.contains(P::inf())); + BOOST_CHECK(!b.contains(P::minus_inf())); + + BOOST_CHECK_EQUAL(b.dimension(), 3); + + b.inflate(2); + auto& bottom = b.get_lower_corner(); + auto& top = b.get_upper_corner(); + BOOST_CHECK_EQUAL(bottom[0], -2); + BOOST_CHECK_EQUAL(bottom[1], -1); + BOOST_CHECK_EQUAL(bottom[2], 0); + BOOST_CHECK_EQUAL(top[0], 6); + BOOST_CHECK_EQUAL(top[1], 7); + BOOST_CHECK_EQUAL(top[2], 8); + + top = P::inf(); + + BOOST_CHECK(b.contains({1, 2, 3})); + BOOST_CHECK(b.contains({1, 8, 3})); + BOOST_CHECK(!b.contains(P::nan())); + BOOST_CHECK(b.contains(P::inf())); + BOOST_CHECK(!b.contains(P::minus_inf())); + + BOOST_CHECK_EQUAL(b.dimension(), 3); + + b.inflate(2); + BOOST_CHECK_EQUAL(bottom[0], -4); + BOOST_CHECK_EQUAL(bottom[1], -3); + BOOST_CHECK_EQUAL(bottom[2], -2); + BOOST_CHECK(top.is_plus_inf()); + + bottom = P::minus_inf(); + top = {4,5,6}; + + BOOST_CHECK(b.contains({1, 2, 3})); + BOOST_CHECK(!b.contains({1, 8, 3})); + BOOST_CHECK(!b.contains(P::nan())); + BOOST_CHECK(!b.contains(P::inf())); + BOOST_CHECK(b.contains(P::minus_inf())); + + BOOST_CHECK_EQUAL(b.dimension(), 3); + + b.inflate(2); + BOOST_CHECK(bottom.is_minus_inf()); + BOOST_CHECK_EQUAL(top[0], 6); + BOOST_CHECK_EQUAL(top[1], 7); + BOOST_CHECK_EQUAL(top[2], 8); + + top = P::inf(); + + BOOST_CHECK(b.contains({1, 2, 3})); + BOOST_CHECK(b.contains({1, 8, 3})); + BOOST_CHECK(!b.contains(P::nan())); + BOOST_CHECK(b.contains(P::inf())); + BOOST_CHECK(b.contains(P::minus_inf())); + + BOOST_CHECK_EQUAL(b.dimension(), 0); + + b.inflate(2); + BOOST_CHECK(bottom.is_minus_inf()); + BOOST_CHECK(top.is_plus_inf()); +} + + + diff --git a/src/Multi_persistence/test/multipersistence_line_unit_test.cpp b/src/Multi_persistence/test/multipersistence_line_unit_test.cpp new file mode 100644 index 0000000000..d86db82d1d --- /dev/null +++ b/src/Multi_persistence/test/multipersistence_line_unit_test.cpp @@ -0,0 +1,154 @@ +/* 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): Hannah Schreiber + * + * Copyright (C) 2024 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#define BOOST_TEST_DYN_LINK +#define BOOST_TEST_MODULE "multi_persistence" +#include +#include + +#include + +using Gudhi::multi_persistence::Line; + +typedef boost::mpl::list list_of_tested_variants; + +BOOST_AUTO_TEST_CASE_TEMPLATE(line_constructors, T, list_of_tested_variants) +{ + using P = typename Line::Point; + + Line l; + BOOST_CHECK(l.base_point().empty()); + BOOST_CHECK(l.direction().empty()); + + P x({1,2,3}); + + Line l2(x); + auto bp = l2.base_point(); + BOOST_CHECK(bp.size() == 3); + BOOST_CHECK_EQUAL(bp[0], 1); + BOOST_CHECK_EQUAL(bp[1], 2); + BOOST_CHECK_EQUAL(bp[2], 3); + BOOST_CHECK(l2.direction().empty()); + + Line l3(std::move(x)); + bp = l3.base_point(); + BOOST_CHECK(bp.size() == 3); + BOOST_CHECK_EQUAL(bp[0], 1); + BOOST_CHECK_EQUAL(bp[1], 2); + BOOST_CHECK_EQUAL(bp[2], 3); + BOOST_CHECK(l3.direction().empty()); + + Line l4({1,2,3}, {4,5,6}); + bp = l4.base_point(); + auto dir = l4.direction(); + BOOST_CHECK(bp.size() == 3); + BOOST_CHECK_EQUAL(bp[0], 1); + BOOST_CHECK_EQUAL(bp[1], 2); + BOOST_CHECK_EQUAL(bp[2], 3); + BOOST_CHECK(dir.size() == 3); + BOOST_CHECK_EQUAL(dir[0], 4); + BOOST_CHECK_EQUAL(dir[1], 5); + BOOST_CHECK_EQUAL(dir[2], 6); + + BOOST_CHECK_THROW(Line l5({1,2,3}, {4,-5,6}), std::invalid_argument); + BOOST_CHECK_THROW(Line l6({1,2,3}, {0,0,0}), std::invalid_argument); + BOOST_CHECK_THROW(Line l7({1,2,3}, {1,2,3,4}), std::invalid_argument); +} + +BOOST_AUTO_TEST_CASE_TEMPLATE(line_intersections, T, list_of_tested_variants) +{ + using P = typename Line::Point; + using KP = typename Line::K_critical_point; + + const double inf = std::numeric_limits::infinity(); + + Line l({1,2,3}, {4,0,6}); + + double t = l.template compute_forward_intersection(P{2,1,3}); + BOOST_CHECK_EQUAL(t, 0.25); + t = l.template compute_forward_intersection(P{2,3,3}); + BOOST_CHECK_EQUAL(t, inf); + + t = l.template compute_forward_intersection(KP({P{2,1,3}, P{2,3,3}})); + BOOST_CHECK_EQUAL(t, 0.25); + + t = l.template compute_backward_intersection(P{2,3,3}); + BOOST_CHECK_EQUAL(t, 0); + t = l.template compute_backward_intersection(P{2,1,3}); + BOOST_CHECK_EQUAL(t, -inf); + + t = l.template compute_backward_intersection(KP({P{2,1,3}, P{2,3,3}})); + BOOST_CHECK_EQUAL(t, 0); + + std::pair bounds = l.get_bounds({{-10, 0, 10}, {10, 4, 10}}); + auto bottom = l[bounds.first]; + auto top = l[bounds.second]; + BOOST_CHECK_EQUAL(bottom.size(), 3); + BOOST_CHECK_EQUAL(bottom[0], T(5. + 2. / 3.)); + BOOST_CHECK_EQUAL(bottom[1], 2); + BOOST_CHECK_EQUAL(bottom[2], T(3. + T(7. / 6.) * 6.)); + BOOST_CHECK_EQUAL(top.size(), 3); + BOOST_CHECK_EQUAL(top[0], T(5. + 2. / 3.)); + BOOST_CHECK_EQUAL(top[1], 2); + BOOST_CHECK_EQUAL(top[2], T(3. + T(7. / 6.) * 6.)); + + bounds = l.get_bounds({{-10, 0, 10}, {10, 1, 10}}); + BOOST_CHECK_EQUAL(bounds.first, P::T_inf); + BOOST_CHECK_EQUAL(bounds.second, -P::T_inf); + BOOST_CHECK(l[bounds.first].is_plus_inf()); + BOOST_CHECK(l[bounds.second].is_minus_inf()); +} + +BOOST_AUTO_TEST_CASE_TEMPLATE(line_other, T, list_of_tested_variants) +{ + using P = typename Line::Point; + + Line l({1,2,3}, {4,0,6}); + + P x = l[1.25]; + BOOST_CHECK_EQUAL(x.size(), 3); + BOOST_CHECK_EQUAL(x[0], 1. + T(1.25) * 4.); + BOOST_CHECK_EQUAL(x[1], 2); + BOOST_CHECK_EQUAL(x[2], 3. + T(1.25) * 6.); + + x = l[0]; + BOOST_CHECK_EQUAL(x.size(), 3); + BOOST_CHECK_EQUAL(x[0], 1); + BOOST_CHECK_EQUAL(x[1], 2); + BOOST_CHECK_EQUAL(x[2], 3); + + x = l[-3.25]; + BOOST_CHECK_EQUAL(x.size(), 3); + BOOST_CHECK_EQUAL(x[0], 1. + T(-3.25) * 4.); + BOOST_CHECK_EQUAL(x[1], 2); + BOOST_CHECK_EQUAL(x[2], 3. + T(-3.25) * 6.); + + l += {2, 5, 6}; + x = l[1.25]; + BOOST_CHECK_EQUAL(x.size(), 3); + BOOST_CHECK_EQUAL(x[0], 3. + T(1.25) * 4.); + BOOST_CHECK_EQUAL(x[1], 7); + BOOST_CHECK_EQUAL(x[2], 9. + T(1.25) * 6.); + + x = l[0]; + BOOST_CHECK_EQUAL(x.size(), 3); + BOOST_CHECK_EQUAL(x[0], 3); + BOOST_CHECK_EQUAL(x[1], 7); + BOOST_CHECK_EQUAL(x[2], 9); + + x = l[-3.25]; + BOOST_CHECK_EQUAL(x.size(), 3); + BOOST_CHECK_EQUAL(x[0], 3. + T(-3.25) * 4.); + BOOST_CHECK_EQUAL(x[1], 7); + BOOST_CHECK_EQUAL(x[2], 9. + T(-3.25) * 6.); +} + + + diff --git a/src/Persistent_cohomology/test/betti_numbers_unit_test.cpp b/src/Persistent_cohomology/test/betti_numbers_unit_test.cpp index f1372a6dc5..c8fddcc4fe 100644 --- a/src/Persistent_cohomology/test/betti_numbers_unit_test.cpp +++ b/src/Persistent_cohomology/test/betti_numbers_unit_test.cpp @@ -19,6 +19,8 @@ struct MiniSTOptions : Gudhi::Simplex_tree_options_minimal { static const bool store_key = true; // I have few vertices typedef short Vertex_handle; + + static const bool is_multi_parameter = false; }; using Mini_simplex_tree = Gudhi::Simplex_tree; diff --git a/src/Persistent_cohomology/test/persistent_cohomology_unit_test.cpp b/src/Persistent_cohomology/test/persistent_cohomology_unit_test.cpp index 5a68ffb600..71c6717087 100644 --- a/src/Persistent_cohomology/test/persistent_cohomology_unit_test.cpp +++ b/src/Persistent_cohomology/test/persistent_cohomology_unit_test.cpp @@ -177,6 +177,7 @@ struct MiniSTOptions { static const bool contiguous_vertices = false; static const bool link_nodes_by_label = false; static const bool stable_simplex_handles = false; + static const bool is_multi_parameter = false; }; using Mini_simplex_tree = Gudhi::Simplex_tree; diff --git a/src/Simplex_tree/benchmark/simplex_tree_cofaces_benchmark.cpp b/src/Simplex_tree/benchmark/simplex_tree_cofaces_benchmark.cpp index d20ea35f2d..56f21ce776 100644 --- a/src/Simplex_tree/benchmark/simplex_tree_cofaces_benchmark.cpp +++ b/src/Simplex_tree/benchmark/simplex_tree_cofaces_benchmark.cpp @@ -91,6 +91,7 @@ struct Stree_basic_cofaces_options { static const bool contiguous_vertices = false; static const bool link_nodes_by_label = false; static const bool stable_simplex_handles = false; + static const bool is_multi_parameter = false; }; struct Stree_fast_cofaces_options : Stree_basic_cofaces_options { @@ -100,6 +101,7 @@ struct Stree_fast_cofaces_options : Stree_basic_cofaces_options { struct Stree_fast_cofaces_stable_simplex_handles_options : Stree_basic_cofaces_options { static const bool link_nodes_by_label = true; static const bool stable_simplex_handles = true; + static const bool is_multi_parameter = false; }; int main(int argc, char *argv[]) { diff --git a/src/Simplex_tree/concept/FiltrationValue.h b/src/Simplex_tree/concept/FiltrationValue.h index 6cf314fa9d..44ff6bc3b2 100644 --- a/src/Simplex_tree/concept/FiltrationValue.h +++ b/src/Simplex_tree/concept/FiltrationValue.h @@ -17,7 +17,19 @@ * its subsimplices of same filtration value) provides an indexing scheme * (see IndexingTag). */ - struct FiltrationValue { - /** \brief Operator < is a StrictWeakOrdering. */ - bool operator<(FiltrationValue f1, FiltrationValue f2); - }; +struct FiltrationValue { + /** \brief Operator < is a StrictWeakOrdering. */ + bool operator<(FiltrationValue f1, FiltrationValue f2); + + /** \brief For multiparameter filtrations, this methods pushes a filtration value + * to the first moment, for operator< such that f < this. For instance, for a one critical filtration, with + * - this = (1,2) + * - x = (2,1) + * + * after calling this method, x should be equal to (2,2). + * This function is called when using, e.g. `make_filtration_non_decreasing`, as the filtration of a simplex + * has to be greater than the filtration of any of its faces. + * */ + void push_to_least_common_upper_bound(const FiltrationValue f); + +}; diff --git a/src/Simplex_tree/concept/SimplexTreeOptions.h b/src/Simplex_tree/concept/SimplexTreeOptions.h index 70f936a2f6..18ad08a3e6 100644 --- a/src/Simplex_tree/concept/SimplexTreeOptions.h +++ b/src/Simplex_tree/concept/SimplexTreeOptions.h @@ -39,5 +39,12 @@ struct SimplexTreeOptions { static const bool link_nodes_by_label; /** @brief If true, Simplex_handle will not be invalidated after insertions or removals. */ static const bool stable_simplex_handles; + /// If true, assumes that SimplexTreeOptions::Filtration_value is vector-like instead of float-like. + /// In that case only, this also assumes that SimplexTreeOptions::Filtration_value is a class, + /// which has a `push_to_least_common_upper_bound` method that allows to push the filtration value `this` onto the set of points + /// \f$ \{ y\in \mathrm{Filtration_value} : y\geq x\}\f$ + /// that are greater than another filtration value \f$ x \f$. + /// An example of such a class is Gudhi::multi_persistence::Finitely_critical_multi_filtration . + static const bool is_multi_parameter; }; diff --git a/src/Simplex_tree/example/CMakeLists.txt b/src/Simplex_tree/example/CMakeLists.txt index 5ea1bbb793..ec24e30c6c 100644 --- a/src/Simplex_tree/example/CMakeLists.txt +++ b/src/Simplex_tree/example/CMakeLists.txt @@ -10,6 +10,10 @@ add_test(NAME Simplex_tree_example_simple_simplex_tree COMMAND $) +add_executable_with_targets ( Simplex_tree_multi_example simplex_tree_multi.cpp TBB::tbb) +add_test(NAME Simplex_tree_multi_example COMMAND $) + + # An example with Simplex-tree using CGAL alpha_shapes_3 if(GMP_FOUND AND TARGET CGAL::CGAL) add_executable ( Simplex_tree_example_alpha_shapes_3_from_off example_alpha_shapes_3_simplex_tree_from_off_file.cpp ) diff --git a/src/Simplex_tree/example/simplex_tree_multi.cpp b/src/Simplex_tree/example/simplex_tree_multi.cpp new file mode 100644 index 0000000000..aa81150cd2 --- /dev/null +++ b/src/Simplex_tree/example/simplex_tree_multi.cpp @@ -0,0 +1,57 @@ +/* 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): David Loiseaux + * + * Copyright (C) 2023 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#include +#include + +#include +#include + +struct ST_MULTI { +public: + typedef Gudhi::linear_indexing_tag Indexing_tag; + typedef int Vertex_handle; + typedef float value_type; + using Filtration_value = Gudhi::multi_filtration::One_critical_filtration; + typedef std::uint32_t Simplex_key; + static const bool store_key = true; + static const bool store_filtration = true; + static const bool contiguous_vertices = false; + static const bool link_nodes_by_label = true; + static const bool stable_simplex_handles = false; + static const bool is_multi_parameter = true; +}; + +using ST = Gudhi::Simplex_tree; + + +int main() { + ST st; + + /* Complex to build. */ + /* 1 */ + /* o */ + /* /X\ */ + /* o---o---o */ + /* 2 0 3 */ + + auto triangle012 = {0, 1, 2}; + auto edge03 = {0, 3}; + st.insert_simplex_and_subfaces(triangle012, {1,2,3}); // {1,2,3} can be any array-like vector-like + st.insert_simplex_and_subfaces(edge03, {4,5,6}); + + auto edge02 = {0, 2}; + ST::Simplex_handle e = st.find(edge02); + // Finitely_critical_multi_filtration has an operator<< + std::cout << st.filtration(e) << std::endl; + GUDHI_CHECK(st.filtration(st.find(edge03)) == ST_MULTI::Filtration_value({4, 5, 6}), + "edge03 does not have the right value."); + +} diff --git a/src/Simplex_tree/include/gudhi/Simplex_tree.h b/src/Simplex_tree/include/gudhi/Simplex_tree.h index 794a1662e5..0d0ee1bf8f 100644 --- a/src/Simplex_tree/include/gudhi/Simplex_tree.h +++ b/src/Simplex_tree/include/gudhi/Simplex_tree.h @@ -162,10 +162,11 @@ class Simplex_tree { Key_simplex_base; struct Filtration_simplex_base_real { - Filtration_simplex_base_real() : filt_(0) {} - Filtration_simplex_base_real(Filtration_value f) : filt_(f) {} - void assign_filtration(Filtration_value f) { filt_ = f; } - Filtration_value filtration() const { return filt_; } + Filtration_simplex_base_real() : filt_{} {} + Filtration_simplex_base_real(const Filtration_value& f) : filt_(f) {} + void assign_filtration(const Filtration_value& f) { filt_ = f; } + const Filtration_value& filtration() const { return filt_; } + Filtration_value& filtration() { return filt_; } private: Filtration_value filt_; }; @@ -177,7 +178,8 @@ class Simplex_tree { void assign_filtration(Filtration_value GUDHI_CHECK_code(f)) { GUDHI_CHECK(f == 0, "filtration value specified for a complex that does not store them"); } - Filtration_value filtration() const { return 0; } + const Filtration_value& filtration() const { return null_value; } + static constexpr Filtration_value null_value={}; }; typedef typename std::conditional::type Filtration_simplex_base; @@ -418,7 +420,9 @@ class Simplex_tree { : null_vertex_(-1), root_(nullptr, null_vertex_), filtration_vect_(), - dimension_(-1) { } + dimension_(-1) { + if constexpr (Options::is_multi_parameter) number_of_parameters_ = 2; + } /** * @brief Construct the simplex tree as the copy of a given simplex tree with eventually different template @@ -457,7 +461,7 @@ class Simplex_tree { * stored in the simplices. * \exception std::invalid_argument In debug mode, if the complex_source is invalid. */ - Simplex_tree(Simplex_tree && complex_source) { + Simplex_tree(Simplex_tree && complex_source) : number_of_parameters_(std::move(complex_source.number_of_parameters_)) { #ifdef DEBUG_TRACES std::clog << "Simplex_tree move constructor" << std::endl; #endif // DEBUG_TRACES @@ -540,6 +544,7 @@ class Simplex_tree { } rec_copy( &root_, &root_source, [](const Filtration_value& fil) -> const Filtration_value& { return fil; }); + if constexpr (Options::is_multi_parameter) number_of_parameters_ = complex_source.number_of_parameters_; } // Copy from complex_source to "this" @@ -563,6 +568,7 @@ class Simplex_tree { } rec_copy(&root_, &root_source, translate_filtration_value); + if constexpr (Options::is_multi_parameter) number_of_parameters_ = complex_source.number_of_parameters_; } /** \brief depth first search, inserts simplices when reaching a leaf. */ @@ -689,7 +695,7 @@ class Simplex_tree { * * Same as `filtration()`, but does not handle `null_simplex()`. */ - static Filtration_value filtration_(Simplex_handle sh) { + static const Filtration_value& filtration_(Simplex_handle sh) { GUDHI_CHECK (sh != null_simplex(), "null simplex"); return sh->second.filtration(); } @@ -724,18 +730,25 @@ class Simplex_tree { * Called on the null_simplex, it returns infinity. * If SimplexTreeOptions::store_filtration is false, returns 0. */ - static Filtration_value filtration(Simplex_handle sh) { + static const Filtration_value& filtration(Simplex_handle sh){ if (sh != null_simplex()) { return sh->second.filtration(); } else { - return std::numeric_limits::infinity(); + return inf_; + } + } + Filtration_value& filtration_mutable(Simplex_handle sh){ + if (sh != null_simplex()) { + return _to_node_it(sh)->second.filtration(); + } else { + return inf_; } } /** \brief Sets the filtration value of a simplex. * \exception std::invalid_argument In debug mode, if sh is a null_simplex. */ - void assign_filtration(Simplex_handle sh, Filtration_value fv) { + void assign_filtration(Simplex_handle sh, const Filtration_value& fv) { GUDHI_CHECK(sh != null_simplex(), std::invalid_argument("Simplex_tree::assign_filtration - cannot assign filtration on null_simplex")); _to_node_it(sh)->second.assign_filtration(fv); @@ -970,14 +983,16 @@ class Simplex_tree { * to the new simplex. * If the insertion fails (the simplex is already there), the bool is set to false. If the insertion * fails and the simplex already in the complex has a filtration value strictly bigger than 'filtration', + * and the simplex tree is not multi parameter (`SimplexTreeOptions::is_multi_parameter == false`), * we assign this simplex with the new value 'filtration', and set the Simplex_handle field of the - * output pair to the Simplex_handle of the simplex. Otherwise, we set the Simplex_handle part to - * null_simplex. - * + * output pair to the Simplex_handle of the simplex. When the simplex tree is multi parameter, + * the existing filtration values are not updated. If the insertion fails for other reasons, + * we set the Simplex_handle part to `null_simplex`. + * */ template > std::pair insert_simplex_raw(const RandomVertexHandleRange& simplex, - Filtration_value filtration) { + const Filtration_value& filtration) { Siblings * curr_sib = &root_; std::pair res_insert; auto vi = simplex.begin(); @@ -1043,7 +1058,7 @@ class Simplex_tree { * .end() return input iterators, with 'value_type' Vertex_handle. */ template> std::pair insert_simplex(const InputVertexRange & simplex, - Filtration_value filtration = 0) { + const Filtration_value& filtration = {}) { auto first = std::begin(simplex); auto last = std::end(simplex); @@ -1072,7 +1087,7 @@ class Simplex_tree { */ template> std::pair insert_simplex_and_subfaces(const InputVertexRange& Nsimplex, - Filtration_value filtration = 0) { + const Filtration_value& filtration = {}) { auto first = std::begin(Nsimplex); auto last = std::end(Nsimplex); @@ -1101,7 +1116,7 @@ class Simplex_tree { std::pair rec_insert_simplex_and_subfaces_sorted(Siblings* sib, ForwardVertexIterator first, ForwardVertexIterator last, - Filtration_value filt) { + const Filtration_value& filt) { // An alternative strategy would be: // - try to find the complete simplex, if found (and low filtration) exit // - insert all the vertices at once in sib @@ -1117,11 +1132,14 @@ class Simplex_tree { // update extra data structures if the insertion is successful update_simplex_tree_after_node_insertion(insertion_result.first); } else { - if (filtration(simplex_one) > filt) { - assign_filtration(simplex_one, filt); - } else { - // FIXME: this interface makes no sense, and it doesn't seem to be tested. - result.first = null_simplex(); + if constexpr (!SimplexTreeOptions::is_multi_parameter) { // Ignores the assign part for multiparameter + // filtrations. + if (filtration(simplex_one) > filt) { + assign_filtration(simplex_one, filt); + } else { + // FIXME: this interface makes no sense, and it doesn't seem to be tested. + result.first = null_simplex(); + } } } if (++first == last) return result; @@ -1484,7 +1502,7 @@ class Simplex_tree { * The complex does not need to be empty before calling this function. However, if a vertex is * already present, its filtration value is not modified, unlike with other insertion functions. */ template - void insert_batch_vertices(VertexRange const& vertices, Filtration_value filt = 0) { + void insert_batch_vertices(VertexRange const& vertices, const Filtration_value& filt ={}) { auto verts = vertices | boost::adaptors::transformed([&](auto v){ return Dit_value_t(v, Node(&root_, filt)); }); root_.members_.insert(boost::begin(verts), boost::end(verts)); @@ -1838,7 +1856,7 @@ class Simplex_tree { static void intersection(std::vector >& intersection, Dictionary_const_it begin1, Dictionary_const_it end1, Dictionary_const_it begin2, Dictionary_const_it end2, - Filtration_value filtration_) { + const Filtration_value& filtration_) { if (begin1 == end1 || begin2 == end2) return; // ----->> while (true) { @@ -2054,18 +2072,33 @@ class Simplex_tree { if (dim == 0) return; // Find the maximum filtration value in the border Boundary_simplex_range&& boundary = boundary_simplex_range(sh); - Boundary_simplex_iterator max_border = std::max_element(std::begin(boundary), std::end(boundary), - [this](Simplex_handle sh1, Simplex_handle sh2) { - return filtration(sh1) < filtration(sh2); - }); + Filtration_value max_filt_border_value; + if constexpr (SimplexTreeOptions::is_multi_parameter) { + // in that case, we assume that Filtration_value has a `push_to_least_common_upper_bound` member to handle this. + max_filt_border_value = Filtration_value(number_of_parameters_); + for (auto& face_sh : boundary) { + max_filt_border_value.push_to_least_common_upper_bound( + filtration(face_sh)); // pushes the value of max_filt_border_value to reach simplex' filtration + } + } else { + Boundary_simplex_iterator max_border = + std::max_element(std::begin(boundary), std::end(boundary), + [this](Simplex_handle sh1, Simplex_handle sh2) { return filtration(sh1) < filtration(sh2); }); + max_filt_border_value = filtration(*max_border); + } - Filtration_value max_filt_border_value = filtration(*max_border); // Replacing if(f=max)) would mean that if f is NaN, we replace it with the max of the children. // That seems more useful than keeping NaN. if (!(sh->second.filtration() >= max_filt_border_value)) { // Store the filtration modification information modified = true; - assign_filtration(sh, max_filt_border_value); + if constexpr (Options::is_multi_parameter){ + auto& to_increase_filtration = filtration_mutable(sh); + to_increase_filtration.push_to_least_common_upper_bound(max_filt_border_value); + } + else{ + assign_filtration(sh, max_filt_border_value); + } } }; // Loop must be from the end to the beginning, as higher dimension simplex are always on the left part of the tree @@ -2094,7 +2127,7 @@ class Simplex_tree { * than it was before. However, `upper_bound_dimension()` will return the old value, which remains a valid upper * bound. If you care, you can call `dimension()` to recompute the exact dimension. */ - bool prune_above_filtration(Filtration_value filtration) { + bool prune_above_filtration(const Filtration_value& filtration) { if (std::numeric_limits::has_infinity && filtration == std::numeric_limits::infinity()) return false; // ---->> bool modified = rec_prune_above_filtration(root(), filtration); @@ -2104,7 +2137,7 @@ class Simplex_tree { } private: - bool rec_prune_above_filtration(Siblings* sib, Filtration_value filt) { + bool rec_prune_above_filtration(Siblings* sib, const Filtration_value& filt) { auto&& list = sib->members(); bool modified = false; bool emptied = false; @@ -2545,7 +2578,7 @@ class Simplex_tree { * @param[in] filt_value The new filtration value. * @param[in] min_dim The minimal dimension. Default value is 0. */ - void reset_filtration(Filtration_value filt_value, int min_dim = 0) { + void reset_filtration(const Filtration_value& filt_value, int min_dim = 0) { rec_reset_filtration(&root_, filt_value, min_dim); clear_filtration(); // Drop the cache. } @@ -2556,7 +2589,7 @@ class Simplex_tree { * @param[in] filt_value The new filtration value. * @param[in] min_depth The minimal depth. */ - void rec_reset_filtration(Siblings * sib, Filtration_value filt_value, int min_depth) { + void rec_reset_filtration(Siblings * sib, const Filtration_value& filt_value, int min_depth) { for (auto sh = sib->members().begin(); sh != sib->members().end(); ++sh) { if (min_depth <= 0) { sh->second.assign_filtration(filt_value); @@ -2728,6 +2761,35 @@ class Simplex_tree { /** \brief Upper bound on the dimension of the simplicial complex.*/ mutable int dimension_; mutable bool dimension_to_be_lowered_ = false; + + // MULTIPERS STUFF + public: + /** + * \brief Sets the number of parameters of the filtrations if SimplexTreeOptions::is_multi_parameter. + * */ + void set_number_of_parameters(int num) { + static_assert(SimplexTreeOptions::is_multi_parameter, + "Cannot set number of parameters of 1-parameter simplextree." + ); + number_of_parameters_ = num; + } + /** + * \brief Gets the number of parameters of the filtrations if SimplexTreeOptions::is_multi_parameter. + * */ + int get_number_of_parameters() const { + if constexpr (SimplexTreeOptions::is_multi_parameter) + return number_of_parameters_; + else + return 1; + } + + // cannot be const due to `filtration_mutable`, TODO : find a proper way to make that const + inline static Filtration_value inf_ = std::numeric_limits::has_infinity ? + std::numeric_limits::infinity() + : std::numeric_limits::max(); /**< Default infinite value. */ + + private: + int number_of_parameters_; /**< Number of parameters of the multi-filtrations when SimplexTreeOptions::is_multi_parameter.-*/ }; // Print a Simplex_tree in os. diff --git a/src/Simplex_tree/include/gudhi/Simplex_tree/Simplex_tree_node_explicit_storage.h b/src/Simplex_tree/include/gudhi/Simplex_tree/Simplex_tree_node_explicit_storage.h index 05bb12b92e..119631f911 100644 --- a/src/Simplex_tree/include/gudhi/Simplex_tree/Simplex_tree_node_explicit_storage.h +++ b/src/Simplex_tree/include/gudhi/Simplex_tree/Simplex_tree_node_explicit_storage.h @@ -43,12 +43,12 @@ struct GUDHI_EMPTY_BASE_CLASS_OPTIMIZATION Simplex_tree_node_explicit_storage typedef typename SimplexTree::Simplex_key Simplex_key; typedef typename SimplexTree::Simplex_data Simplex_data; - Simplex_tree_node_explicit_storage(Siblings* sib = nullptr, Filtration_value filtration = 0) + Simplex_tree_node_explicit_storage(Siblings* sib = nullptr, const Filtration_value& filtration = {}) : SimplexTree::Filtration_simplex_base(filtration), children_(sib) {} //will fail to compile when called with SimplexTree::Options::store_key is false. - Simplex_tree_node_explicit_storage(Siblings* sib, Filtration_value filtration, Simplex_key key) + Simplex_tree_node_explicit_storage(Siblings* sib, const Filtration_value& filtration, Simplex_key key) : SimplexTree::Filtration_simplex_base(filtration), SimplexTree::Key_simplex_base(key), children_(sib) {} 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 d8db37d244..5a38b75bc4 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 @@ -39,6 +39,7 @@ struct Simplex_tree_options_default { static const bool contiguous_vertices = false; static const bool link_nodes_by_label = false; static const bool stable_simplex_handles = false; + static const bool is_multi_parameter = false; }; /** Model of SimplexTreeOptions. @@ -55,6 +56,7 @@ struct Simplex_tree_options_full_featured { static const bool contiguous_vertices = false; static const bool link_nodes_by_label = true; static const bool stable_simplex_handles = true; + static const bool is_multi_parameter = false; }; /** Model of SimplexTreeOptions. @@ -71,6 +73,7 @@ struct Simplex_tree_options_minimal { static const bool contiguous_vertices = false; static const bool link_nodes_by_label = false; static const bool stable_simplex_handles = false; + static const bool is_multi_parameter = false; }; /** @private @brief Model of SimplexTreeOptions, faster than `Simplex_tree_options_default` but note the unsafe @@ -88,6 +91,7 @@ struct Simplex_tree_options_fast_persistence { static const bool contiguous_vertices = true; static const bool link_nodes_by_label = false; static const bool stable_simplex_handles = false; + static const bool is_multi_parameter = false; }; /** @}*/ // end addtogroup simplex_tree diff --git a/src/Simplex_tree/test/CMakeLists.txt b/src/Simplex_tree/test/CMakeLists.txt index 6a1ecac7f6..df51d7d2cb 100644 --- a/src/Simplex_tree/test/CMakeLists.txt +++ b/src/Simplex_tree/test/CMakeLists.txt @@ -6,6 +6,9 @@ file(COPY "simplex_tree_for_unit_test.txt" DESTINATION ${CMAKE_CURRENT_BINARY_DI add_executable_with_targets(Simplex_tree_test_unit simplex_tree_unit_test.cpp TBB::tbb) gudhi_add_boost_test(Simplex_tree_test_unit) +add_executable_with_targets(Simplex_tree_multi_test_unit simplex_tree_multi_unit_test.cpp TBB::tbb) +gudhi_add_boost_test(Simplex_tree_multi_test_unit) + add_executable_with_targets(Simplex_tree_remove_test_unit simplex_tree_remove_unit_test.cpp TBB::tbb) gudhi_add_boost_test(Simplex_tree_remove_test_unit) diff --git a/src/Simplex_tree/test/simplex_tree_ctor_and_move_unit_test.cpp b/src/Simplex_tree/test/simplex_tree_ctor_and_move_unit_test.cpp index a5a7f5e1e6..077419adf7 100644 --- a/src/Simplex_tree/test/simplex_tree_ctor_and_move_unit_test.cpp +++ b/src/Simplex_tree/test/simplex_tree_ctor_and_move_unit_test.cpp @@ -187,6 +187,7 @@ struct Simplex_tree_options_custom_fil_values_default { static const bool contiguous_vertices = false; static const bool link_nodes_by_label = false; static const bool stable_simplex_handles = false; + static const bool is_multi_parameter = false; }; struct Simplex_tree_options_custom_fil_values_fast_persistence { @@ -199,6 +200,7 @@ struct Simplex_tree_options_custom_fil_values_fast_persistence { static const bool contiguous_vertices = true; static const bool link_nodes_by_label = false; static const bool stable_simplex_handles = false; + static const bool is_multi_parameter = false; }; struct Simplex_tree_options_custom_fil_values_full_featured { @@ -211,6 +213,7 @@ struct Simplex_tree_options_custom_fil_values_full_featured { static const bool contiguous_vertices = false; static const bool link_nodes_by_label = true; static const bool stable_simplex_handles = true; + static const bool is_multi_parameter = false; }; typedef boost::mpl::list, diff --git a/src/Simplex_tree/test/simplex_tree_edge_expansion_unit_test.cpp b/src/Simplex_tree/test/simplex_tree_edge_expansion_unit_test.cpp index 618bde1b34..f655629686 100644 --- a/src/Simplex_tree/test/simplex_tree_edge_expansion_unit_test.cpp +++ b/src/Simplex_tree/test/simplex_tree_edge_expansion_unit_test.cpp @@ -35,6 +35,7 @@ struct Simplex_tree_options_fast_cofaces { static const bool contiguous_vertices = false; static const bool link_nodes_by_label = true; static const bool stable_simplex_handles = false; + static const bool is_multi_parameter = false; }; struct Simplex_tree_options_stable_simplex_handles_fast_cofaces : Simplex_tree_options_fast_cofaces { diff --git a/src/Simplex_tree/test/simplex_tree_graph_expansion_unit_test.cpp b/src/Simplex_tree/test/simplex_tree_graph_expansion_unit_test.cpp index e523cce1ec..e9c4584690 100644 --- a/src/Simplex_tree/test/simplex_tree_graph_expansion_unit_test.cpp +++ b/src/Simplex_tree/test/simplex_tree_graph_expansion_unit_test.cpp @@ -31,6 +31,7 @@ struct Simplex_tree_options_stable_simplex_handles { static const bool contiguous_vertices = false; static const bool link_nodes_by_label = false; static const bool stable_simplex_handles = true; + static const bool is_multi_parameter = false; }; typedef boost::mpl::list, diff --git a/src/Simplex_tree/test/simplex_tree_iostream_operator_unit_test.cpp b/src/Simplex_tree/test/simplex_tree_iostream_operator_unit_test.cpp index edeb40ca50..8cb3c11798 100644 --- a/src/Simplex_tree/test/simplex_tree_iostream_operator_unit_test.cpp +++ b/src/Simplex_tree/test/simplex_tree_iostream_operator_unit_test.cpp @@ -25,6 +25,7 @@ struct MyOptions : Simplex_tree_options_minimal { // Not doing persistence, so we don't need those static const bool store_key = false; static const bool store_filtration = false; + static const bool is_multi_parameter = false; // I have few vertices typedef short Vertex_handle; }; diff --git a/src/Simplex_tree/test/simplex_tree_make_filtration_non_decreasing_unit_test.cpp b/src/Simplex_tree/test/simplex_tree_make_filtration_non_decreasing_unit_test.cpp index 34adee3b98..7faf4e94a9 100644 --- a/src/Simplex_tree/test/simplex_tree_make_filtration_non_decreasing_unit_test.cpp +++ b/src/Simplex_tree/test/simplex_tree_make_filtration_non_decreasing_unit_test.cpp @@ -33,6 +33,7 @@ struct Simplex_tree_options_stable_simplex_handles { static const bool contiguous_vertices = false; static const bool link_nodes_by_label = false; static const bool stable_simplex_handles = true; + static const bool is_multi_parameter = false; }; typedef boost::mpl::list, diff --git a/src/Simplex_tree/test/simplex_tree_multi_unit_test.cpp b/src/Simplex_tree/test/simplex_tree_multi_unit_test.cpp new file mode 100644 index 0000000000..09d21ca145 --- /dev/null +++ b/src/Simplex_tree/test/simplex_tree_multi_unit_test.cpp @@ -0,0 +1,982 @@ +/* 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): David Loiseaux, Vincent Rouvreau + * + * Copyright (C) 2023 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#include +#include +#include +#include // std::pair, std::make_pair +#include // float comparison +#include +#include // greater +#include // std::tie +#include // for std::distance +#include // for std::size_t + +#define BOOST_TEST_DYN_LINK +#define BOOST_TEST_MODULE "simplex_tree_multi" +#include +#include + +#include +#include +#include + +using namespace Gudhi; +using OneCriticalFiltration = Gudhi::multi_filtration::One_critical_filtration; +using Stree_options_multi_filtration = Gudhi::multi_persistence::Simplex_tree_options_multidimensional_filtration; + +using typeST_STD = Simplex_tree; +using Stree_multi = Simplex_tree; + +typedef boost::mpl::list list_of_tested_variants; + +template +void test_empty_simplex_tree(typeST& tst) { + typedef typename typeST::Vertex_handle Vertex_handle; + const Vertex_handle DEFAULT_VERTEX_VALUE = Vertex_handle(- 1); + BOOST_CHECK(tst.null_vertex() == DEFAULT_VERTEX_VALUE); + BOOST_CHECK(tst.num_vertices() == (size_t) 0); + BOOST_CHECK(tst.num_simplices() == (size_t) 0); + BOOST_CHECK(tst.is_empty()); + BOOST_CHECK(tst.num_simplices_by_dimension() == std::vector()); + typename typeST::Siblings* STRoot = tst.root(); + BOOST_CHECK(STRoot != nullptr); + BOOST_CHECK(STRoot->oncles() == nullptr); + BOOST_CHECK(STRoot->parent() == DEFAULT_VERTEX_VALUE); + BOOST_CHECK(tst.dimension() == -1); +} + +template +void test_iterators_on_empty_simplex_tree(typeST& tst) { + std::clog << "Iterator on vertices: " << std::endl; + for (auto vertex : tst.complex_vertex_range()) { + std::clog << "vertice:" << vertex << std::endl; + BOOST_CHECK(false); // shall be empty + } + std::clog << "Iterator on simplices: " << std::endl; + for (auto simplex : tst.complex_simplex_range()) { + BOOST_CHECK(simplex != simplex); // shall be empty - to remove warning of non-used simplex + } + + std::clog + << "Iterator on Simplices in the filtration, with [filtration value]:" + << std::endl; + for (auto f_simplex : tst.filtration_simplex_range()) { + BOOST_CHECK(false); // shall be empty + std::clog << "test_iterators_on_empty_simplex_tree - filtration=" + << tst.filtration(f_simplex) << std::endl; + } +} + +BOOST_AUTO_TEST_CASE_TEMPLATE(simplex_tree_when_empty, typeST, list_of_tested_variants) { + typedef std::pair typePairSimplexBool; + typedef std::vector typeVectorVertex; + + std::clog << "********************************************************************" << std::endl; + std::clog << "TEST OF DEFAULT CONSTRUCTOR" << std::endl; + typeST st; + + test_empty_simplex_tree(st); + + test_iterators_on_empty_simplex_tree(st); + // TEST OF EMPTY INSERTION + std::clog << "TEST OF EMPTY INSERTION" << std::endl; + typeVectorVertex simplexVectorEmpty; + BOOST_CHECK(simplexVectorEmpty.empty() == true); + typePairSimplexBool returnEmptyValue = st.insert_simplex(simplexVectorEmpty, 0.0); + BOOST_CHECK(returnEmptyValue.first == typeST::null_simplex()); + BOOST_CHECK(returnEmptyValue.second == true); + + test_empty_simplex_tree(st); + + test_iterators_on_empty_simplex_tree(st); +} + +bool AreAlmostTheSame(Gudhi::multi_filtration::One_critical_filtration a, + Gudhi::multi_filtration::One_critical_filtration b) { + assert(a.size() == b.size()); + for (auto i=0u; i std::numeric_limits::epsilon()) + return false; + return true; +} + + +template +void test_simplex_tree_contains(typeST& simplexTree, typeSimplex& simplex, int pos) { + auto f_simplex = simplexTree.filtration_simplex_range().begin() + pos; + + std::clog << "test_simplex_tree_contains - filtration=" << simplexTree.filtration(*f_simplex) << "||" << simplex.second << std::endl; + BOOST_CHECK(AreAlmostTheSame(simplexTree.filtration(*f_simplex), simplex.second)); + + int simplexIndex = simplex.first.size() - 1; + std::sort(simplex.first.begin(), simplex.first.end()); // if the simplex wasn't sorted, the next test could fail + for (auto vertex : simplexTree.simplex_vertex_range(*f_simplex)) { + std::clog << "test_simplex_tree_contains - vertex=" << vertex << "||" << simplex.first.at(simplexIndex) << std::endl; + BOOST_CHECK(vertex == simplex.first.at(simplexIndex)); + BOOST_CHECK(simplexIndex >= 0); + simplexIndex--; + } +} + +template +void test_simplex_tree_insert_returns_true(const typePairSimplexBool& returnValue) { + BOOST_CHECK(returnValue.second == true); + // Simplex_handle = boost::container::flat_map< typeST::Vertex_handle, Node >::iterator + typename typeST::Simplex_handle shReturned = returnValue.first; + BOOST_CHECK(shReturned != typeST::null_simplex()); +} + +// Global variables +int dim_max = -1; + +template +void set_and_test_simplex_tree_dim_fil(typeST& simplexTree, int vectorSize, const Filtration_value& fil) { + if (vectorSize > dim_max + 1) { + dim_max = vectorSize - 1; + simplexTree.set_dimension(dim_max); + std::clog << " set_and_test_simplex_tree_dim_fil - dim_max=" << dim_max + << std::endl; + } + + BOOST_CHECK(simplexTree.dimension() == dim_max); + + // Another way to count simplices: + size_t num_simp = 0; + for (auto f_simplex : simplexTree.complex_simplex_range()) { + // Remove warning + (void) f_simplex; + num_simp++; + } + + BOOST_CHECK(simplexTree.num_simplices() == num_simp); +} + +BOOST_AUTO_TEST_CASE_TEMPLATE(simplex_tree_insertion, typeST, list_of_tested_variants) { + typedef typename typeST::Filtration_value Filtration_value; + typedef std::pair typePairSimplexBool; + typedef std::vector typeVectorVertex; + typedef std::pair typeSimplex; + const Filtration_value FIRST_FILTRATION_VALUE = {0.1,0.15}; + const Filtration_value SECOND_FILTRATION_VALUE = {0.2, 0.25}; + const Filtration_value THIRD_FILTRATION_VALUE = {0.3, 0.35}; + const Filtration_value FOURTH_FILTRATION_VALUE = {0.4, 0.45}; + // reset since we run the test several times + dim_max = -1; + + // TEST OF INSERTION + std::clog << "********************************************************************" << std::endl; + std::clog << "TEST OF INSERTION" << std::endl; + typeST st; + + // ++ FIRST + std::clog << " - INSERT 0" << std::endl; + typeVectorVertex firstSimplexVector{0}; + BOOST_CHECK(firstSimplexVector.size() == 1); + typeSimplex firstSimplex = std::make_pair(firstSimplexVector, Filtration_value(FIRST_FILTRATION_VALUE)); + typePairSimplexBool returnValue = st.insert_simplex(firstSimplex.first, firstSimplex.second); + + test_simplex_tree_insert_returns_true(returnValue); + set_and_test_simplex_tree_dim_fil(st, firstSimplexVector.size(), firstSimplex.second); + BOOST_CHECK(st.num_vertices() == (size_t) 1); + + // ++ SECOND + std::clog << " - INSERT 1" << std::endl; + typeVectorVertex secondSimplexVector{1}; + BOOST_CHECK(secondSimplexVector.size() == 1); + typeSimplex secondSimplex = std::make_pair(secondSimplexVector, Filtration_value(FIRST_FILTRATION_VALUE)); + returnValue = st.insert_simplex(secondSimplex.first, secondSimplex.second); + + test_simplex_tree_insert_returns_true(returnValue); + set_and_test_simplex_tree_dim_fil(st, secondSimplexVector.size(), secondSimplex.second); + BOOST_CHECK(st.num_vertices() == (size_t) 2); + + // ++ THIRD + std::clog << " - INSERT (0,1)" << std::endl; + typeVectorVertex thirdSimplexVector{0, 1}; + BOOST_CHECK(thirdSimplexVector.size() == 2); + typeSimplex thirdSimplex = std::make_pair(thirdSimplexVector, Filtration_value(SECOND_FILTRATION_VALUE)); + returnValue = st.insert_simplex(thirdSimplex.first, thirdSimplex.second); + + test_simplex_tree_insert_returns_true(returnValue); + set_and_test_simplex_tree_dim_fil(st, thirdSimplexVector.size(), thirdSimplex.second); + BOOST_CHECK(st.num_vertices() == (size_t) 2); // Not incremented !! + + // ++ FOURTH + std::clog << " - INSERT 2" << std::endl; + typeVectorVertex fourthSimplexVector{2}; + BOOST_CHECK(fourthSimplexVector.size() == 1); + typeSimplex fourthSimplex = std::make_pair(fourthSimplexVector, Filtration_value(FIRST_FILTRATION_VALUE)); + returnValue = st.insert_simplex(fourthSimplex.first, fourthSimplex.second); + + test_simplex_tree_insert_returns_true(returnValue); + set_and_test_simplex_tree_dim_fil(st, fourthSimplexVector.size(), fourthSimplex.second); + BOOST_CHECK(st.num_vertices() == (size_t) 3); + + // ++ FIFTH + std::clog << " - INSERT (2,0)" << std::endl; + typeVectorVertex fifthSimplexVector{2, 0}; + BOOST_CHECK(fifthSimplexVector.size() == 2); + typeSimplex fifthSimplex = std::make_pair(fifthSimplexVector, Filtration_value(SECOND_FILTRATION_VALUE)); + returnValue = st.insert_simplex(fifthSimplex.first, fifthSimplex.second); + + test_simplex_tree_insert_returns_true(returnValue); + set_and_test_simplex_tree_dim_fil(st, fifthSimplexVector.size(), fifthSimplex.second); + BOOST_CHECK(st.num_vertices() == (size_t) 3); // Not incremented !! + + // ++ SIXTH + std::clog << " - INSERT (2,1)" << std::endl; + typeVectorVertex sixthSimplexVector{2, 1}; + BOOST_CHECK(sixthSimplexVector.size() == 2); + typeSimplex sixthSimplex = std::make_pair(sixthSimplexVector, Filtration_value(SECOND_FILTRATION_VALUE)); + returnValue = st.insert_simplex(sixthSimplex.first, sixthSimplex.second); + + test_simplex_tree_insert_returns_true(returnValue); + set_and_test_simplex_tree_dim_fil(st, sixthSimplexVector.size(), sixthSimplex.second); + BOOST_CHECK(st.num_vertices() == (size_t) 3); // Not incremented !! + + // ++ SEVENTH + std::clog << " - INSERT (2,1,0)" << std::endl; + typeVectorVertex seventhSimplexVector{2, 1, 0}; + BOOST_CHECK(seventhSimplexVector.size() == 3); + typeSimplex seventhSimplex = std::make_pair(seventhSimplexVector, Filtration_value(THIRD_FILTRATION_VALUE)); + returnValue = st.insert_simplex(seventhSimplex.first, seventhSimplex.second); + + test_simplex_tree_insert_returns_true(returnValue); + set_and_test_simplex_tree_dim_fil(st, seventhSimplexVector.size(), seventhSimplex.second); + BOOST_CHECK(st.num_vertices() == (size_t) 3); // Not incremented !! + + // ++ EIGHTH + std::clog << " - INSERT 3" << std::endl; + typeVectorVertex eighthSimplexVector{3}; + BOOST_CHECK(eighthSimplexVector.size() == 1); + typeSimplex eighthSimplex = std::make_pair(eighthSimplexVector, Filtration_value(FIRST_FILTRATION_VALUE)); + returnValue = st.insert_simplex(eighthSimplex.first, eighthSimplex.second); + + test_simplex_tree_insert_returns_true(returnValue); + set_and_test_simplex_tree_dim_fil(st, eighthSimplexVector.size(), eighthSimplex.second); + BOOST_CHECK(st.num_vertices() == (size_t) 4); + + // ++ NINTH + std::clog << " - INSERT (3,0)" << std::endl; + typeVectorVertex ninethSimplexVector{3, 0}; + BOOST_CHECK(ninethSimplexVector.size() == 2); + typeSimplex ninethSimplex = std::make_pair(ninethSimplexVector, Filtration_value(SECOND_FILTRATION_VALUE)); + returnValue = st.insert_simplex(ninethSimplex.first, ninethSimplex.second); + + test_simplex_tree_insert_returns_true(returnValue); + set_and_test_simplex_tree_dim_fil(st, ninethSimplexVector.size(), ninethSimplex.second); + BOOST_CHECK(st.num_vertices() == (size_t) 4); // Not incremented !! + + // ++ TENTH + std::clog << " - INSERT 0 (already inserted)" << std::endl; + typeVectorVertex tenthSimplexVector{0}; + BOOST_CHECK(tenthSimplexVector.size() == 1); + // With a different filtration value + typeSimplex tenthSimplex = std::make_pair(tenthSimplexVector, Filtration_value(FOURTH_FILTRATION_VALUE)); + returnValue = st.insert_simplex(tenthSimplex.first, tenthSimplex.second); + + BOOST_CHECK(returnValue.second == false); + // Simplex_handle = boost::container::flat_map< typeST::Vertex_handle, Node >::iterator + typename typeST::Simplex_handle shReturned = returnValue.first; + BOOST_CHECK(shReturned == typeST::null_simplex()); + std::clog << "st.num_vertices()=" << st.num_vertices() << std::endl; + BOOST_CHECK(st.num_vertices() == (size_t) 4); // Not incremented !! + BOOST_CHECK(st.dimension() == dim_max); + + // ++ ELEVENTH + std::clog << " - INSERT (2,1,0) (already inserted)" << std::endl; + typeVectorVertex eleventhSimplexVector{2, 1, 0}; + BOOST_CHECK(eleventhSimplexVector.size() == 3); + typeSimplex eleventhSimplex = std::make_pair(eleventhSimplexVector, Filtration_value(FOURTH_FILTRATION_VALUE)); + returnValue = st.insert_simplex(eleventhSimplex.first, eleventhSimplex.second); + + BOOST_CHECK(returnValue.second == false); + // Simplex_handle = boost::container::flat_map< typeST::Vertex_handle, Node >::iterator + shReturned = returnValue.first; + BOOST_CHECK(shReturned == typeST::null_simplex()); + BOOST_CHECK(st.num_vertices() == (size_t) 4); // Not incremented !! + BOOST_CHECK(st.dimension() == dim_max); + BOOST_CHECK(st.num_simplices_by_dimension() == std::vector({4, 4, 1})); + + /* Inserted simplex: */ + /* 1 */ + /* o */ + /* /X\ */ + /* o---o---o */ + /* 2 0 3 */ + + // [0.1] 0 + // [0.1] 1 + // [0.1] 2 + // [0.1] 3 + // [0.2] 1 0 + // [0.2] 2 0 + // [0.2] 2 1 + // [0.2] 3 0 + // [0.3] 2 1 0 + // !! Be careful, simplex are sorted by filtration value on insertion !! + std::clog << "simplex_tree_insertion - first - 0" << std::endl; + test_simplex_tree_contains(st, firstSimplex, 0); // (0) -> 0 + std::clog << "simplex_tree_insertion - second - 1" << std::endl; + test_simplex_tree_contains(st, secondSimplex, 1); // (1) -> 1 + std::clog << "simplex_tree_insertion - third - 4" << std::endl; + test_simplex_tree_contains(st, thirdSimplex, 4); // (0,1) -> 4 + std::clog << "simplex_tree_insertion - fourth - 2" << std::endl; + test_simplex_tree_contains(st, fourthSimplex, 2); // (2) -> 2 + std::clog << "simplex_tree_insertion - fifth - 5" << std::endl; + test_simplex_tree_contains(st, fifthSimplex, 5); // (2,0) -> 5 + std::clog << "simplex_tree_insertion - sixth - 6" << std::endl; + test_simplex_tree_contains(st, sixthSimplex, 6); //(2,1) -> 6 + std::clog << "simplex_tree_insertion - seventh - 8" << std::endl; + test_simplex_tree_contains(st, seventhSimplex, 8); // (2,1,0) -> 8 + std::clog << "simplex_tree_insertion - eighth - 3" << std::endl; + test_simplex_tree_contains(st, eighthSimplex, 3); // (3) -> 3 + std::clog << "simplex_tree_insertion - ninth - 7" << std::endl; + test_simplex_tree_contains(st, ninethSimplex, 7); // (3,0) -> 7 + + // Display the Simplex_tree - Can not be done in the middle of 2 inserts + std::clog << "The complex contains " << st.num_simplices() << " simplices" << std::endl; + std::clog << " - dimension " << st.dimension() << std::endl; + std::clog << std::endl << std::endl << "Iterator on Simplices in the filtration, with [filtration value]:" << std::endl; + for (auto f_simplex : st.filtration_simplex_range()) { + std::clog << " " << "[" << st.filtration(f_simplex) << "] "; + for (auto vertex : st.simplex_vertex_range(f_simplex)) { + std::clog << (int) vertex << " "; + } + std::clog << std::endl; + } + +} + +BOOST_AUTO_TEST_CASE_TEMPLATE(NSimplexAndSubfaces_tree_insertion, typeST, list_of_tested_variants) { + typedef std::pair typePairSimplexBool; + typedef std::vector typeVectorVertex; + typedef typename typeST::Filtration_value F; + typedef std::pair typeSimplex; + std::clog << "********************************************************************" << std::endl; + std::clog << "TEST OF RECURSIVE INSERTION" << std::endl; + typeST st; + typePairSimplexBool returnValue; + int position = 0; + + // ++ FIRST + std::clog << " - INSERT (2,1,0)" << std::endl; + typeVectorVertex SimplexVector1{2, 1, 0}; + BOOST_CHECK(SimplexVector1.size() == 3); + returnValue = st.insert_simplex_and_subfaces(SimplexVector1); + + BOOST_CHECK(st.num_vertices() == (size_t) 3); // +3 (2, 1 and 0 are not existing) + + // Check it is well inserted + BOOST_CHECK(true == returnValue.second); + position = 0; + std::sort(SimplexVector1.begin(), SimplexVector1.end(), std::greater()); + for (auto vertex : st.simplex_vertex_range(returnValue.first)) { + // Check returned Simplex_handle + std::clog << "vertex = " << vertex << " | vector[" << position << "] = " << SimplexVector1[position] << std::endl; + BOOST_CHECK(vertex == SimplexVector1[position]); + position++; + } + + // ++ SECOND + std::clog << " - INSERT 3" << std::endl; + typeVectorVertex SimplexVector2{3}; + BOOST_CHECK(SimplexVector2.size() == 1); + returnValue = st.insert_simplex_and_subfaces(SimplexVector2); + + BOOST_CHECK(st.num_vertices() == (size_t) 4); // +1 (3 is not existing) + + // Check it is well inserted + BOOST_CHECK(true == returnValue.second); + position = 0; + std::sort(SimplexVector2.begin(), SimplexVector2.end(), std::greater()); + for (auto vertex : st.simplex_vertex_range(returnValue.first)) { + // Check returned Simplex_handle + std::clog << "vertex = " << vertex << " | vector[" << position << "] = " << SimplexVector2[position] << std::endl; + BOOST_CHECK(vertex == SimplexVector2[position]); + position++; + } + + // ++ THIRD + std::clog << " - INSERT (0,3)" << std::endl; + typeVectorVertex SimplexVector3{3, 0}; + BOOST_CHECK(SimplexVector3.size() == 2); + returnValue = st.insert_simplex_and_subfaces(SimplexVector3); + + BOOST_CHECK(st.num_vertices() == (size_t) 4); // Not incremented (all are existing) + + // Check it is well inserted + BOOST_CHECK(true == returnValue.second); + position = 0; + std::sort(SimplexVector3.begin(), SimplexVector3.end(), std::greater()); + for (auto vertex : st.simplex_vertex_range(returnValue.first)) { + // Check returned Simplex_handle + std::clog << "vertex = " << vertex << " | vector[" << position << "] = " << SimplexVector3[position] << std::endl; + BOOST_CHECK(vertex == SimplexVector3[position]); + position++; + } + + // ++ FOURTH + std::clog << " - INSERT (1,0) (already inserted)" << std::endl; + typeVectorVertex SimplexVector4{1, 0}; + BOOST_CHECK(SimplexVector4.size() == 2); + returnValue = st.insert_simplex_and_subfaces(SimplexVector4); + + BOOST_CHECK(st.num_vertices() == (size_t) 4); // Not incremented (all are existing) + + // Check it was not inserted (already there from {2,1,0} insertion) + BOOST_CHECK(false == returnValue.second); + + // ++ FIFTH + std::clog << " - INSERT (3,4,5)" << std::endl; + typeVectorVertex SimplexVector5{3, 4, 5}; + BOOST_CHECK(SimplexVector5.size() == 3); + returnValue = st.insert_simplex_and_subfaces(SimplexVector5); + + BOOST_CHECK(st.num_vertices() == (size_t) 6); + + // Check it is well inserted + BOOST_CHECK(true == returnValue.second); + position = 0; + std::sort(SimplexVector5.begin(), SimplexVector5.end(), std::greater()); + for (auto vertex : st.simplex_vertex_range(returnValue.first)) { + // Check returned Simplex_handle + std::clog << "vertex = " << vertex << " | vector[" << position << "] = " << SimplexVector5[position] << std::endl; + BOOST_CHECK(vertex == SimplexVector5[position]); + position++; + } + + // ++ SIXTH + std::clog << " - INSERT (0,1,6,7)" << std::endl; + typeVectorVertex SimplexVector6{0, 1, 6, 7}; + BOOST_CHECK(SimplexVector6.size() == 4); + returnValue = st.insert_simplex_and_subfaces(SimplexVector6); + + BOOST_CHECK(st.num_vertices() == (size_t) 8); // +2 (6 and 7 are not existing - 0 and 1 are already existing) + + // Check it is well inserted + BOOST_CHECK(true == returnValue.second); + position = 0; + std::sort(SimplexVector6.begin(), SimplexVector6.end(), std::greater()); + for (auto vertex : st.simplex_vertex_range(returnValue.first)) { + // Check returned Simplex_handle + std::clog << "vertex = " << vertex << " | vector[" << position << "] = " << SimplexVector6[position] << std::endl; + BOOST_CHECK(vertex == SimplexVector6[position]); + position++; + } + + /* Inserted simplex: */ + /* 1 6 */ + /* o---o */ + /* /X\7/ */ + /* o---o---o---o */ + /* 2 0 3\X/4 */ + /* o */ + /* 5 */ + /* */ + /* In other words: */ + /* A facet [2,1,0] */ + /* An edge [0,3] */ + /* A facet [3,4,5] */ + /* A cell [0,1,6,7] */ + + typeSimplex simplexPair1 = std::make_pair(SimplexVector1, F()); + typeSimplex simplexPair2 = std::make_pair(SimplexVector2, F()); + typeSimplex simplexPair3 = std::make_pair(SimplexVector3, F()); + typeSimplex simplexPair4 = std::make_pair(SimplexVector4, F()); + typeSimplex simplexPair5 = std::make_pair(SimplexVector5, F()); + typeSimplex simplexPair6 = std::make_pair(SimplexVector6, F()); + test_simplex_tree_contains(st, simplexPair1, 6); // (2,1,0) is in position 6 + test_simplex_tree_contains(st, simplexPair2, 7); // (3) is in position 7 + test_simplex_tree_contains(st, simplexPair3, 8); // (3,0) is in position 8 + test_simplex_tree_contains(st, simplexPair4, 2); // (1,0) is in position 2 + test_simplex_tree_contains(st, simplexPair5, 14); // (3,4,5) is in position 14 + test_simplex_tree_contains(st, simplexPair6, 26); // (7,6,1,0) is in position 26 + + // ------------------------------------------------------------------------------------------------------------------ + // Find in the simplex_tree + // ------------------------------------------------------------------------------------------------------------------ + typeVectorVertex simpleSimplexVector{1}; + typename typeST::Simplex_handle simplexFound = st.find(simpleSimplexVector); + std::clog << "**************IS THE SIMPLEX {1} IN THE SIMPLEX TREE ?\n"; + if (simplexFound != st.null_simplex()) + std::clog << "***+ YES IT IS!\n"; + else + std::clog << "***- NO IT ISN'T\n"; + // Check it is found + BOOST_CHECK(simplexFound != st.null_simplex()); + + typeVectorVertex unknownSimplexVector{15}; + simplexFound = st.find(unknownSimplexVector); + std::clog << "**************IS THE SIMPLEX {15} IN THE SIMPLEX TREE ?\n"; + if (simplexFound != st.null_simplex()) + std::clog << "***+ YES IT IS!\n"; + else + std::clog << "***- NO IT ISN'T\n"; + // Check it is NOT found + BOOST_CHECK(simplexFound == st.null_simplex()); + + simplexFound = st.find(SimplexVector6); + std::clog << "**************IS THE SIMPLEX {0,1,6,7} IN THE SIMPLEX TREE ?\n"; + if (simplexFound != st.null_simplex()) + std::clog << "***+ YES IT IS!\n"; + else + std::clog << "***- NO IT ISN'T\n"; + // Check it is found + BOOST_CHECK(simplexFound != st.null_simplex()); + + typeVectorVertex otherSimplexVector{1, 15}; + simplexFound = st.find(otherSimplexVector); + std::clog << "**************IS THE SIMPLEX {15,1} IN THE SIMPLEX TREE ?\n"; + if (simplexFound != st.null_simplex()) + std::clog << "***+ YES IT IS!\n"; + else + std::clog << "***- NO IT ISN'T\n"; + // Check it is NOT found + BOOST_CHECK(simplexFound == st.null_simplex()); + + typeVectorVertex invSimplexVector{1, 2, 0}; + simplexFound = st.find(invSimplexVector); + std::clog << "**************IS THE SIMPLEX {1,2,0} IN THE SIMPLEX TREE ?\n"; + if (simplexFound != st.null_simplex()) + std::clog << "***+ YES IT IS!\n"; + else + std::clog << "***- NO IT ISN'T\n"; + // Check it is found + BOOST_CHECK(simplexFound != st.null_simplex()); + + // Display the Simplex_tree - Can not be done in the middle of 2 inserts + std::clog << "The complex contains " << st.num_simplices() << " simplices" << std::endl; + std::clog << " - dimension " << st.dimension() << std::endl; + std::clog << std::endl << std::endl << "Iterator on Simplices in the filtration, with [filtration value]:" << std::endl; + for (auto f_simplex : st.filtration_simplex_range()) { + std::clog << " " << "[" << st.filtration(f_simplex) << "] "; + for (auto vertex : st.simplex_vertex_range(f_simplex)) { + std::clog << (int) vertex << " "; + } + std::clog << std::endl; + } +} + + +template +void test_simplex_is_vertex(typeST& st, typename typeST::Simplex_handle sh, typename typeST::Vertex_handle v) { + BOOST_CHECK(st.dimension(sh) == 0); + auto&& r = st.simplex_vertex_range(sh); + auto i = std::begin(r); + BOOST_CHECK(*i == v); + BOOST_CHECK(++i == std::end(r)); +} + +BOOST_AUTO_TEST_CASE_TEMPLATE(simplex_tree_reset_filtration, typeST, list_of_tested_variants) { + std::clog << "********************************************************************" << std::endl; + std::clog << "TEST RESET FILTRATION" << std::endl; + typeST st; + + st.insert_simplex_and_subfaces({2, 1, 0}, OneCriticalFiltration({2.,1.})); + st.insert_simplex_and_subfaces({3, 0}, OneCriticalFiltration({1.,2.})); + st.insert_simplex_and_subfaces({3, 4, 5}, OneCriticalFiltration({3.,4.})); + st.insert_simplex_and_subfaces({0, 1, 6, 7}, OneCriticalFiltration({4.,3.})); + std::cout <<"TRUC "<< st.filtration(st.find({2,1,0})) << std::endl; + /* Inserted simplex: */ + /* 1 6 */ + /* o---o */ + /* /X\7/ */ + /* o---o---o---o */ + /* 2 0 3\X/4 */ + /* o */ + /* 5 */ + + for (auto f_simplex : st.skeleton_simplex_range(3)) { + std::clog << "vertex = ("; + for (auto vertex : st.simplex_vertex_range(f_simplex)) { + std::clog << vertex << ","; + } + std::clog << ") - filtration = " << st.filtration(f_simplex); + std::clog << " - dimension = " << st.dimension(f_simplex) << std::endl; + // Guaranteed by construction + BOOST_CHECK(st.filtration(f_simplex) >= OneCriticalFiltration({1.,1.})); + } + + // dimension until 5 even if simplex tree is of dimension 3 to test the limits + for(int dimension = 5; dimension >= 0; dimension --) { + std::clog << "### reset_filtration - dimension = " << dimension << "\n"; + st.reset_filtration(st.inf_, dimension); + for (auto f_simplex : st.skeleton_simplex_range(3)) { + std::clog << "vertex = ("; + for (auto vertex : st.simplex_vertex_range(f_simplex)) { + std::clog << vertex << ","; + } + std::clog << ") - filtration = " << st.filtration(f_simplex); + std::clog << " - dimension = " << st.dimension(f_simplex) << std::endl; + if (st.dimension(f_simplex) < dimension) + BOOST_CHECK(st.filtration(f_simplex) >= OneCriticalFiltration({1.,1})); + else + BOOST_CHECK(st.filtration(f_simplex) == st.inf_); + } + } + +} + +BOOST_AUTO_TEST_CASE_TEMPLATE(simplex_tree_clear, typeST, list_of_tested_variants) { + std::clog << "********************************************************************" << std::endl; + std::clog << "TEST SIMPLEX TREE CLEAR" << std::endl; + typeST st; + st.insert_simplex_and_subfaces({0, 1}, OneCriticalFiltration({1.5})); + st.initialize_filtration(); + st.clear(); + BOOST_CHECK(st.num_vertices() == 0); + BOOST_CHECK(st.num_simplices() == 0); + BOOST_CHECK(st.upper_bound_dimension() == -1); + BOOST_CHECK(st.dimension() == -1); + BOOST_CHECK(boost::size(st.filtration_simplex_range()) == 0); + typeST st_empty; + BOOST_CHECK(st == st_empty); + st.insert_simplex_and_subfaces({0}, OneCriticalFiltration({2.5})); + BOOST_CHECK(boost::size(st.cofaces_simplex_range(st.find({0}), 1)) == 0); +} + +BOOST_AUTO_TEST_CASE_TEMPLATE(multify_simplex_tree, typeST, list_of_tested_variants) { + std::clog << "********************************************************************" << std::endl; + std::clog << "TEST MULTIFY FLATTEN LINEAR PROJECTION" << std::endl; + + typeST_STD st; + typeST st_multi; + st.insert_simplex_and_subfaces({1,2,3}, 1.); + BOOST_CHECK(st.get_number_of_parameters() == 1); + int num_parameters = 3; + Gudhi::multi_persistence::multify(st,st_multi,num_parameters,{2.,3.}); //fills st_multi by simplices of filtration [{st.filtration(sh)}, default_values[0],default_values[1], ...] + BOOST_CHECK(st_multi.get_number_of_parameters() == num_parameters); // num parameters is defined by multify + BOOST_CHECK(st_multi.num_simplices() == st.num_simplices()); // simplicial complexes should be the same + for (auto sh : st_multi.complex_simplex_range()){ + const auto& filtration = st_multi.filtration(sh); + BOOST_CHECK(filtration == OneCriticalFiltration({1,2,3})); // Checks the filtration values + } + std::clog << "********************************************************************" << std::endl; + std::clog << "TEST FLATTEN" << std::endl; + + for (int dimension=0; dimension < 3; dimension++){ + st.clear(); + Gudhi::multi_persistence::flatten(st, st_multi, dimension); + for (auto sh : st.complex_simplex_range()){ + BOOST_CHECK(st.filtration(sh) == dimension +1); + } + } + + // std::clog << "********************************************************************" << std::endl; + // std::clog << "TEST LINEAR PROJECTION" << std::endl; + // // st has already the same simplicial complex as st_multi, and the filtration values of st_multi are all {1,2,3} + // Gudhi::multi_persistence::linear_projection(st,st_multi,{17,37,73}); // sets the filtration values of st to the dot product of st_multi.filtration and {17,37,73}. + // for (auto sh : st.complex_simplex_range()){ + // BOOST_CHECK(st.filtration(sh) == 1*17 + 2*37 + 3*73); + // } +} + + + + +BOOST_AUTO_TEST_CASE(simplex_tree_multi_assign_filtration) { + std::clog << "********************************************************************" << std::endl; + std::clog << "TEST MULTI FILTRATION INSERT SIMPLEX AND SUBFACES" << std::endl; + Stree_multi st; + + typename Stree_multi::Simplex_handle sh; + bool success = false; + const OneCriticalFiltration multi_filt_1 = {1., 2.}; + std::tie(sh, success) = st.insert_simplex_and_subfaces({0, 1}, multi_filt_1); + BOOST_CHECK(success); + BOOST_CHECK(sh != st.null_simplex()); + // Only [0,1], [0] and [1] are already inserted + const OneCriticalFiltration multi_filt_2 = {3., 2., 1.}; + std::tie(sh, success) = st.insert_simplex_and_subfaces({2, 1, 0}, multi_filt_2); + BOOST_CHECK(success); + BOOST_CHECK(sh != st.null_simplex()); + // Already inserted + std::tie(sh, success) = st.insert_simplex_and_subfaces({0, 2}, {4.}); + BOOST_CHECK(!success); + BOOST_CHECK(sh != st.null_simplex()); + + // Check filtration values of an already inserted simplex + sh = st.find({1, 2}); + BOOST_CHECK(sh != st.null_simplex()); + BOOST_CHECK(st.filtration(sh) == multi_filt_2); + // And assign a new value + OneCriticalFiltration const multi_filt_3 = {5.}; + st.assign_filtration(sh, multi_filt_3); + + for (auto f_simplex : st.complex_simplex_range()) { + std::clog << "vertex = ("; + for (auto vertex : st.simplex_vertex_range(f_simplex)) { + std::clog << vertex << ","; + } + std::clog << ") - filtration = " << st.filtration(f_simplex); + std::clog << " - dimension = " << st.dimension(f_simplex) << std::endl; + } + + // Check all filtration values + sh = st.find({2, 1, 0}); + BOOST_CHECK(sh != st.null_simplex()); + BOOST_CHECK(st.filtration(sh) == multi_filt_2); + sh = st.find({1, 0}); + BOOST_CHECK(sh != st.null_simplex()); + BOOST_CHECK(st.filtration(sh) == multi_filt_1); + sh = st.find({2, 0}); + BOOST_CHECK(sh != st.null_simplex()); + BOOST_CHECK(st.filtration(sh) == multi_filt_2); + sh = st.find({0}); + BOOST_CHECK(sh != st.null_simplex()); + BOOST_CHECK(st.filtration(sh) == multi_filt_1); + sh = st.find({2, 1}); + BOOST_CHECK(sh != st.null_simplex()); + BOOST_CHECK(st.filtration(sh) == multi_filt_3); + sh = st.find({1}); + BOOST_CHECK(sh != st.null_simplex()); + BOOST_CHECK(st.filtration(sh) == multi_filt_1); + sh = st.find({2}); + BOOST_CHECK(sh != st.null_simplex()); + BOOST_CHECK(st.filtration(sh) == multi_filt_2); + +} + +BOOST_AUTO_TEST_CASE(simplex_tree_multi_reset_filtration) { + std::clog << "********************************************************************" << std::endl; + std::clog << "TEST RESET MULTI FILTRATION" << std::endl; + Stree_multi st; + + st.insert_simplex_and_subfaces({2, 1, 0}, {3., 2.}); + st.insert_simplex_and_subfaces({3, 0}, {2., 3.}); + st.insert_simplex_and_subfaces({3, 4, 5}, {3., 2.}); + st.insert_simplex_and_subfaces({0, 1, 6, 7}, {4.,4.}); + + /* Inserted simplex: */ + /* 1 6 */ + /* o---o */ + /* /X\7/ */ + /* o---o---o---o */ + /* 2 0 3\X/4 */ + /* o */ + /* 5 */ + + for (auto f_simplex : st.complex_simplex_range()) { + std::clog << "vertex = ("; + for (auto vertex : st.simplex_vertex_range(f_simplex)) { + std::clog << vertex << ","; + } + std::clog << ") - filtration = " << st.filtration(f_simplex); + std::clog << " - dimension = " << st.dimension(f_simplex) << std::endl; + // Guaranteed by construction + BOOST_CHECK(st.filtration(f_simplex) >= 2.); + } + + OneCriticalFiltration new_filt = {0., 1., 2.}; + // dimension until 5 even if simplex tree is of dimension 3 to test the limits + for(int dimension = 5; dimension >= 0; dimension --) { + std::clog << "### reset_filtration - dimension = " << dimension << "\n"; + st.reset_filtration(new_filt, dimension); + for (auto f_simplex : st.complex_simplex_range()) { + std::clog << "vertex = ("; + for (auto vertex : st.simplex_vertex_range(f_simplex)) { + std::clog << vertex << ","; + } + std::clog << ") - filtration = " << st.filtration(f_simplex); + std::clog << " - dimension = " << st.dimension(f_simplex) << std::endl; + if (st.dimension(f_simplex) >= dimension) + BOOST_CHECK(st.filtration(f_simplex) == new_filt); + } + } +} + +BOOST_AUTO_TEST_CASE(simplex_tree_multi_filtration_multify) { + std::clog << "********************************************************************" << std::endl; + std::clog << "TEST MULTI FILTRATION MULTIFY" << std::endl; + + Simplex_tree<> st; + st.insert_simplex_and_subfaces({2, 1, 0}, 3.); + st.insert_simplex_and_subfaces({3, 0}, 2.); + st.insert_simplex_and_subfaces({3, 4, 5}, 3.); + st.insert_simplex_and_subfaces({0, 1, 6, 7}, 4.); + + /* Inserted simplex: */ + /* 1 6 */ + /* o---o */ + /* /X\7/ */ + /* o---o---o---o */ + /* 2 0 3\X/4 */ + /* o */ + /* 5 */ + + Stree_multi st_multi; + OneCriticalFiltration default_multi (3, std::numeric_limits::quiet_NaN()); + Gudhi::multi_persistence::multify(st, st_multi, 3, default_multi); + + for (auto f_simplex : st_multi.complex_simplex_range()) { + std::clog << "vertex = ("; + for (auto vertex : st_multi.simplex_vertex_range(f_simplex)) { + std::clog << vertex << ","; + } + auto multi_filtration = st_multi.filtration(f_simplex); + std::clog << ") - filtration = " << multi_filtration << std::endl; + BOOST_CHECK(!std::isnan(multi_filtration[0])); + BOOST_CHECK(std::isnan(multi_filtration[1])); + BOOST_CHECK(std::isnan(multi_filtration[2])); + } + + { + Simplex_tree<> copy; + Gudhi::multi_persistence::flatten(copy, st_multi, 0); + BOOST_CHECK(st == copy); + } + + { + Simplex_tree<> copy; + Gudhi::multi_persistence::flatten(copy, st_multi, 1); + for (auto f_simplex : copy.complex_simplex_range()) { + std::clog << "vertex = ("; + for (auto vertex : copy.simplex_vertex_range(f_simplex)) { + std::clog << vertex << ","; + } + std::clog << ") - filtration = " << copy.filtration(f_simplex) << std::endl; + BOOST_CHECK(std::isnan(copy.filtration(f_simplex))); + } + } + +} + +BOOST_AUTO_TEST_CASE(simplex_tree_multi_filtration_numeric_limits) { + std::clog << "********************************************************************" << std::endl; + std::clog << "TEST MULTI FILTRATION NUMERIC LIMITS" << std::endl; + + // NaN + auto nan_multi = std::numeric_limits::quiet_NaN(); + BOOST_CHECK(nan_multi.size() == 1); + BOOST_CHECK(std::isnan(nan_multi[0])); + std::clog << nan_multi << std::endl; + + // Inf + auto inf_multi = std::numeric_limits::infinity(); + BOOST_CHECK(inf_multi.size() == 1); + BOOST_CHECK(std::isinf(inf_multi[0])); + std::clog << inf_multi << std::endl; +} + +BOOST_AUTO_TEST_CASE(make_filtration_non_decreasing_on_multi) { + std::clog << "********************************************************************" << std::endl; + std::clog << "TEST MULTI FILTRATION NON DECREASING" << std::endl; + Stree_multi st; + + st.insert_simplex_and_subfaces({2, 1, 0}, {3., 4.}); + st.insert_simplex_and_subfaces({3, 0}, {2., 3.}); + st.insert_simplex_and_subfaces({3, 4, 5}, {3., 4.}); + + /* Inserted simplex: */ + /* 1 */ + /* o */ + /* /X\ */ + /* o---o---o---o */ + /* 2 0 3\X/4 */ + /* o */ + /* 5 */ + + + for (auto f_simplex : st.complex_simplex_range()) { + std::clog << "vertex = ("; + for (auto vertex : st.simplex_vertex_range(f_simplex)) { + std::clog << vertex << ","; + } + std::clog << ") - filtration = " << st.filtration(f_simplex) << std::endl; + } + auto filt = st.filtration(st.find({3, 0})); + std::clog << "filtration([3,0]) = " << filt << std::endl; + BOOST_CHECK(filt == OneCriticalFiltration({2., 3.})); + + st.set_number_of_parameters(2); + st.make_filtration_non_decreasing(); + + filt = st.filtration(st.find({3, 0})); + std::clog << "filtration([3,0]) = " << filt << std::endl; + BOOST_CHECK(filt == OneCriticalFiltration({3., 4.})); + + st.assign_filtration(st.find({3, 0}), OneCriticalFiltration({2.9, 4.})); + filt = st.filtration(st.find({3, 0})); + std::clog << "filtration([3,0]) = " << filt << std::endl; + BOOST_CHECK(filt == OneCriticalFiltration({2.9, 4.})); + + st.make_filtration_non_decreasing(); + + filt = st.filtration(st.find({3, 0})); + std::clog << "filtration([3,0]) = " << filt << std::endl; + BOOST_CHECK(filt == OneCriticalFiltration({3., 4.})); + + st.assign_filtration(st.find({3, 0}), OneCriticalFiltration({5., 3.99})); + filt = st.filtration(st.find({3, 0})); + std::clog << "filtration([3,0]) = " << filt << std::endl; + BOOST_CHECK(filt == OneCriticalFiltration({5., 3.99})); + + st.make_filtration_non_decreasing(); + + filt = st.filtration(st.find({3, 0})); + std::clog << "filtration([3,0]) = " << filt << std::endl; + BOOST_CHECK(filt == OneCriticalFiltration({5., 4.})); +} + +BOOST_AUTO_TEST_CASE(make_filtration_non_decreasing_on_multi_nan_values) { + Stree_multi st; + + BOOST_CHECK(std::numeric_limits::quiet_NaN().is_nan()); + BOOST_CHECK(std::numeric_limits::infinity().is_plus_inf()); + + st.insert_simplex_and_subfaces({2, 1, 0}, {1.,2.,3.}); + st.insert_simplex_and_subfaces({3, 0}, {1.,2.,3.}); + st.insert_simplex_and_subfaces({3, 4, 5}, {1.,2.,3.}); + + st.assign_filtration(st.find({0}), std::numeric_limits::quiet_NaN()); + st.assign_filtration(st.find({3}), std::numeric_limits::infinity()); + + /* Inserted simplex: */ + /* 1 */ + /* o */ + /* /X\ */ + /* o---o---o---o */ + /* 2 0 3\X/4 */ + /* o */ + /* 5 */ + + // Default number of parameter is 2 + BOOST_CHECK(st.get_number_of_parameters() == 2); + + st.set_number_of_parameters(3); + BOOST_CHECK(st.get_number_of_parameters() == 3); + + std::clog << "SPECIFIC CASE:" << std::endl; + std::clog << "Insertion with NaN values does not ensure the filtration values are non decreasing" << std::endl; + st.make_filtration_non_decreasing(); + + std::clog << "Check that NaN filtrations are ignored, and inf filtrations are propagated." << std::endl; + for (auto f_simplex : st.complex_simplex_range()) { + auto filt = st.filtration(f_simplex); + bool contains3 = false; + bool contains0 = false; + int dim = -1; + std::clog << "Simplex "; + for (auto vertex : st.simplex_vertex_range(f_simplex)){ + std::clog << vertex << " "; + if (vertex == 3) contains3 = true; + else if (vertex == 0) contains0 = true; + dim++; + } + bool is_zero = dim == 0 && contains0; + std::clog << "Filtration: " << filt << std::endl; + if (is_zero) BOOST_CHECK(filt.is_nan()); + else if (contains3) BOOST_CHECK(filt.is_plus_inf()); + else BOOST_CHECK(filt == OneCriticalFiltration({1.,2.,3.})); + } +} + diff --git a/src/Simplex_tree/test/simplex_tree_remove_unit_test.cpp b/src/Simplex_tree/test/simplex_tree_remove_unit_test.cpp index f2527050bd..a7a4c49f64 100644 --- a/src/Simplex_tree/test/simplex_tree_remove_unit_test.cpp +++ b/src/Simplex_tree/test/simplex_tree_remove_unit_test.cpp @@ -24,6 +24,7 @@ struct MyOptions : Simplex_tree_options_minimal { // Not doing persistence, so we don't need those static const bool store_key = false; static const bool store_filtration = false; + static const bool is_multi_parameter = false; // I have few vertices typedef short Vertex_handle; }; @@ -567,4 +568,4 @@ BOOST_AUTO_TEST_CASE(mini_prune_above_dimension) { BOOST_CHECK(st.num_simplices() == 0); BOOST_CHECK(st.upper_bound_dimension() == -1); BOOST_CHECK(st.dimension() == -1); -} \ No newline at end of file +} diff --git a/src/Simplex_tree/test/simplex_tree_serialization_unit_test.cpp b/src/Simplex_tree/test/simplex_tree_serialization_unit_test.cpp index a8450cfa72..4259333389 100644 --- a/src/Simplex_tree/test/simplex_tree_serialization_unit_test.cpp +++ b/src/Simplex_tree/test/simplex_tree_serialization_unit_test.cpp @@ -38,6 +38,8 @@ struct Low_options : Gudhi::Simplex_tree_options_full_featured { typedef std::uint8_t Vertex_handle; // Maximum number of simplices to compute persistence is 2^8 - 1 = 255. One is reserved for null_key typedef std::uint8_t Simplex_key; + + static const bool is_multi_parameter = false; }; struct Stable_options : Gudhi::Simplex_tree_options_default { diff --git a/src/cmake/modules/GUDHI_third_party_libraries.cmake b/src/cmake/modules/GUDHI_third_party_libraries.cmake index 250401ff3a..088ab1b8e5 100644 --- a/src/cmake/modules/GUDHI_third_party_libraries.cmake +++ b/src/cmake/modules/GUDHI_third_party_libraries.cmake @@ -200,4 +200,4 @@ if (WITH_GUDHI_PYTHON) option(WITH_GUDHI_PYTHON_RUNTIME_LIBRARY_DIRS "Build with setting runtime_library_dirs. Useful when setting rpath is not allowed" ON) -endif (WITH_GUDHI_PYTHON) +endif (WITH_GUDHI_PYTHON) \ No newline at end of file diff --git a/src/common/doc/examples.h b/src/common/doc/examples.h index d363024beb..33f8683d0d 100644 --- a/src/common/doc/examples.h +++ b/src/common/doc/examples.h @@ -18,6 +18,7 @@ * @example simple_simplex_tree.cpp * @example simplex_tree_from_cliques_of_graph.cpp * @example example_alpha_shapes_3_simplex_tree_from_off_file.cpp + * @example simplex_tree_multi.cpp * \section Persistent_cohomology_example_section Persistent_cohomology * @example custom_persistence_sort.cpp * @example rips_persistence_step_by_step.cpp @@ -141,4 +142,3 @@ * \section Persistence_fields_example_section Persistence_field * @example example_field_operations.cpp */ - diff --git a/src/python/CMakeLists.txt b/src/python/CMakeLists.txt index 4be5ce9602..8314661f04 100644 --- a/src/python/CMakeLists.txt +++ b/src/python/CMakeLists.txt @@ -715,4 +715,4 @@ if(MATPLOTLIB_FOUND) endif() # Set missing or not modules -set(GUDHI_MODULES ${GUDHI_MODULES} "python" CACHE INTERNAL "GUDHI_MODULES") +set(GUDHI_MODULES ${GUDHI_MODULES} "python" CACHE INTERNAL "GUDHI_MODULES") \ No newline at end of file diff --git a/src/python/include/Simplex_tree_interface.h b/src/python/include/Simplex_tree_interface.h index c162feb0cd..cade7aecfe 100644 --- a/src/python/include/Simplex_tree_interface.h +++ b/src/python/include/Simplex_tree_interface.h @@ -38,6 +38,7 @@ struct Simplex_tree_options_for_python { static const bool contiguous_vertices = false; static const bool link_nodes_by_label = false; static const bool stable_simplex_handles = false; + static const bool is_multi_parameter = false; };