From 60b7325bacc7c103b817560e7167ceb8ae115865 Mon Sep 17 00:00:00 2001 From: David Loiseaux Date: Thu, 28 Sep 2023 16:04:13 +0200 Subject: [PATCH] CPP part only --- .../test/betti_numbers_unit_test.cpp | 2 + .../test/persistent_cohomology_unit_test.cpp | 1 + .../simplex_tree_cofaces_benchmark.cpp | 1 + src/Simplex_tree/concept/SimplexTreeOptions.h | 2 + src/Simplex_tree/include/gudhi/Simplex_tree.h | 93 +++++-- .../gudhi/Simplex_tree/Simplex_tree_multi.h | 231 ++++++++++++++++++ .../Simplex_tree_node_explicit_storage.h | 2 +- .../Simplex_tree/multi_filtrations/box.h | 164 +++++++++++++ .../finitely_critical_filtrations.h | 150 ++++++++++++ .../Simplex_tree/multi_filtrations/line.h | 94 +++++++ .../simplex_tree_ctor_and_move_unit_test.cpp | 1 + ...simplex_tree_graph_expansion_unit_test.cpp | 1 + ...ke_filtration_non_decreasing_unit_test.cpp | 1 + .../test/simplex_tree_unit_test.cpp | 1 + 14 files changed, 720 insertions(+), 24 deletions(-) create mode 100644 src/Simplex_tree/include/gudhi/Simplex_tree/Simplex_tree_multi.h create mode 100644 src/Simplex_tree/include/gudhi/Simplex_tree/multi_filtrations/box.h create mode 100644 src/Simplex_tree/include/gudhi/Simplex_tree/multi_filtrations/finitely_critical_filtrations.h create mode 100644 src/Simplex_tree/include/gudhi/Simplex_tree/multi_filtrations/line.h diff --git a/src/Persistent_cohomology/test/betti_numbers_unit_test.cpp b/src/Persistent_cohomology/test/betti_numbers_unit_test.cpp index eccd17b21f..96b880cc94 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_full_featured { 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 f6e001b631..e4353977b6 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 Simplex_tree_options_stable_simplex_handles { 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; }; int main(int argc, char *argv[]) { diff --git a/src/Simplex_tree/concept/SimplexTreeOptions.h b/src/Simplex_tree/concept/SimplexTreeOptions.h index 0c95d3011c..45e5f03bca 100644 --- a/src/Simplex_tree/concept/SimplexTreeOptions.h +++ b/src/Simplex_tree/concept/SimplexTreeOptions.h @@ -31,5 +31,7 @@ struct SimplexTreeOptions { static const bool link_nodes_by_label; /// If true, Simplex_handle will not be invalidated after insertions or removals. static const bool stable_simplex_handles; + /// If true, assumes that Filtration_value is vector-like instead of float-like. This also assumes that Filtration_values is a class, which has a push_to method to push a filtration value $x$ onto $this>=0$. + static const bool is_multi_parameter; }; diff --git a/src/Simplex_tree/include/gudhi/Simplex_tree.h b/src/Simplex_tree/include/gudhi/Simplex_tree.h index 9f93672d63..0effd69a34 100644 --- a/src/Simplex_tree/include/gudhi/Simplex_tree.h +++ b/src/Simplex_tree/include/gudhi/Simplex_tree.h @@ -151,16 +151,18 @@ class Simplex_tree { Key_simplex_base; struct Filtration_simplex_base_real { - Filtration_simplex_base_real() : filt_(0) {} - void assign_filtration(Filtration_value f) { filt_ = f; } - Filtration_value filtration() const { return filt_; } + Filtration_simplex_base_real() : filt_{} {} + 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_; }; struct Filtration_simplex_base_dummy { Filtration_simplex_base_dummy() {} 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; @@ -576,7 +578,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(); } @@ -604,18 +606,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_; + } + } + static Filtration_value& filtration_mutable(Simplex_handle sh){ + if (sh != null_simplex()) { + return 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")); sh->second.assign_filtration(fv); @@ -829,7 +838,7 @@ class Simplex_tree { */ 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(); @@ -895,7 +904,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); @@ -924,7 +933,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); @@ -953,7 +962,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 @@ -970,8 +979,20 @@ class Simplex_tree { Simplex_handle simplex_one = insertion_result.first; bool one_is_new = insertion_result.second; if (!one_is_new) { - if (filtration(simplex_one) > filt) { - assign_filtration(simplex_one, filt); + if (!(filtration(simplex_one) <= filt)) { + if constexpr (SimplexTreeOptions::is_multi_parameter){ + // By default, does nothing and assumes that the user is smart. + if (filt < filtration(simplex_one)){ + // placeholder for comparable filtrations + } + else{ + // placeholder for incomparable filtrations + } + } + else{ // non-multiparameter + assign_filtration(simplex_one, filt); + } + } else { // FIXME: this interface makes no sense, and it doesn't seem to be tested. insertion_result.first = null_simplex(); @@ -1316,7 +1337,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)); @@ -1403,7 +1424,7 @@ class Simplex_tree { static void intersection(std::vector >& intersection, Dictionary_it begin1, Dictionary_it end1, Dictionary_it begin2, Dictionary_it end2, - Filtration_value filtration_) { + const Filtration_value& filtration_) { if (begin1 == end1 || begin2 == end2) return; // ----->> while (true) { @@ -1610,12 +1631,22 @@ 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), + Filtration_value max_filt_border_value; + if constexpr (SimplexTreeOptions::is_multi_parameter){ + // in that case, we assume that Filtration_value has a `push_to` member to handle this. + max_filt_border_value = Filtration_value(this->number_of_parameters_); + for (auto &face_sh : boundary){ + max_filt_border_value.push_to(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), [](Simplex_handle sh1, Simplex_handle sh2) { return filtration(sh1) < filtration(sh2); }); - - Filtration_value max_filt_border_value = filtration(*max_border); + 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)) { @@ -1650,7 +1681,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); @@ -1660,7 +1691,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(); auto last = std::remove_if(list.begin(), list.end(), [this,filt](Dit_value_t& simplex) { if (simplex.second.filtration() <= filt) return false; @@ -2070,7 +2101,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. } @@ -2081,7 +2112,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); @@ -2246,6 +2277,18 @@ class Simplex_tree { /** \brief Upper bound on the dimension of the simplicial complex.*/ int dimension_; bool dimension_to_be_lowered_ = false; + +//MULTIPERS STUFF +public: + void set_number_of_parameters(int num){ + number_of_parameters_ = num; + } + int get_number_of_parameters() const{ + return number_of_parameters_; + } + inline static Filtration_value inf_ = std::numeric_limits::infinity(); +private: + int number_of_parameters_; }; // Print a Simplex_tree in os. @@ -2297,6 +2340,7 @@ struct Simplex_tree_options_full_featured { 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, faster than `Simplex_tree_options_full_featured` but note the unsafe @@ -2314,6 +2358,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; }; /** Model of SimplexTreeOptions, faster cofaces than `Simplex_tree_options_full_featured`, note the @@ -2331,6 +2376,8 @@ 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; + }; /** @}*/ // end addtogroup simplex_tree diff --git a/src/Simplex_tree/include/gudhi/Simplex_tree/Simplex_tree_multi.h b/src/Simplex_tree/include/gudhi/Simplex_tree/Simplex_tree_multi.h new file mode 100644 index 0000000000..b2eae92cab --- /dev/null +++ b/src/Simplex_tree/include/gudhi/Simplex_tree/Simplex_tree_multi.h @@ -0,0 +1,231 @@ +/* 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) 2014 Inria """ + * + * + * Modification(s): + * - 2022/11 David Loiseaux / Hannah Schreiber : added multify / flatten to interface standard simplextree. + * - YYYY/MM Author: Description of the modification + */ +#ifndef SIMPLEX_TREE_MULTI_H_ +#define SIMPLEX_TREE_MULTI_H_ + +#include +#include +#include +#include + + + + + +namespace Gudhi::multiparameter { +/** Model of SimplexTreeOptions. + * + * Maximum number of simplices to compute persistence is std::numeric_limits::max() + * (about 4 billions of simplices). */ + + +struct Simplex_tree_options_multidimensional_filtration { +public: + typedef linear_indexing_tag Indexing_tag; + typedef int Vertex_handle; + typedef float value_type; + using Filtration_value = multi_filtrations::Finitely_critical_multi_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 options_multi = Simplex_tree_options_multidimensional_filtration; +using options_std = Simplex_tree_options_full_featured; +using simplextree_std = Simplex_tree; +using simplextree_multi = Simplex_tree; + +using multi_filtration_type = std::vector; +using multi_filtration_grid = std::vector; + +// Turns a 1-parameter simplextree into a multiparameter simplextree, and keeps the 1-filtration in the 1st axis +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); + for (auto i = 0u; i(default_values.size()), static_cast(num_parameters-1));i++) + 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) + f[0] = st.filtration(simplex_handle); + auto filtration_copy = f; + st_multi.insert_simplex(simplex,filtration_copy); + + } +} + + + +// Turns a multi-parameter simplextree into a 1-parameter simplextree +template +void flatten(simplextree_std &st, simplextree_multi &st_multi, const int dimension = 0){ + for (const auto &simplex_handle : st_multi.complex_simplex_range()){ + std::vector simplex; + for (auto vertex : st_multi.simplex_vertex_range(simplex_handle)) + simplex.push_back(vertex); + typename simplextree_multi::Options::value_type f = dimension >= 0 ? st_multi.filtration(simplex_handle)[dimension] : 0; + st.insert_simplex(simplex,f); + } +} + +// Applies a linear form (i.e. scalar product with Riesz rpz) to the filtration to flatten a simplextree multi +template +void linear_projection(simplextree_std &st, simplextree_multi &st_multi, const std::vector& linear_form){ + auto sh = st.complex_simplex_range().begin(); + auto sh_multi = st_multi.complex_simplex_range().begin(); + auto end = st.complex_simplex_range().end(); + typename simplextree_multi::Options::Filtration_value multi_filtration; + for (; sh != end; ++sh, ++sh_multi){ + multi_filtration = st_multi.filtration(*sh_multi); + auto projected_filtration = multi_filtration.linear_projection(linear_form); + st.assign_filtration(*sh, projected_filtration); + } +} + +template +void flatten_diag(simplextree_std &st, simplextree_multi &st_multi, const std::vector basepoint, int dimension){ + multi_filtrations::Line l(basepoint); + for (const auto &simplex_handle : st_multi.complex_simplex_range()){ + std::vector simplex; + for (auto vertex : st_multi.simplex_vertex_range(simplex_handle)) + simplex.push_back(vertex); + + std::vector f = st_multi.filtration(simplex_handle); + if (dimension <0) dimension = 0; + typename simplextree_multi::Options::value_type new_filtration = l.push_forward(f)[dimension]; + st.insert_simplex(simplex,new_filtration); + } +} + + + + + +/// @brief turns filtration value x into coordinates in the grid +/// @tparam out_type +/// @param x +/// @param grid +/// @return +template +inline void find_coordinates(vector_like& x, const multi_filtration_grid &grid){ + // TODO: optimize with, e.g., dichotomy + + for (auto parameter = 0u; parameter < grid.size(); parameter++){ + const auto& filtration = grid[parameter]; // assumes its sorted + const auto to_project = x[parameter]; + if constexpr (std::numeric_limits::has_infinity) + if (to_project == std::numeric_limits::infinity()){ + x[parameter] = std::numeric_limits::infinity(); + continue; + } + if (to_project >= filtration.back()){ + x[parameter] = filtration.size()-1; + continue; + } // deals with infinite value at the end of the grid + + unsigned int i = 0; + while (to_project > filtration[i] && i &st_multi = *(Gudhi::Simplex_tree*)(splxptr); + auto num_parameters = static_cast(st_multi.get_number_of_parameters()); + if (grid.size() != num_parameters){ + std::cerr << "Bad grid !" << std::endl; + throw; + } + for (const auto &simplex_handle : st_multi.complex_simplex_range()){ + auto& simplex_filtration = st_multi.filtration_mutable(simplex_handle); + find_coordinates(simplex_filtration, grid); // turns the simplexfiltration into coords in the grid + if (!coordinate_values){ + for (auto parameter = 0u; parameter < num_parameters; parameter++) + simplex_filtration[parameter] = grid[parameter][simplex_filtration[parameter]]; + } + } + return; +} + +// retrieves the filtration values of a simplextree. Useful to generate a grid. +std::vector get_filtration_values(const uintptr_t splxptr, const std::vector °rees){ + Simplex_tree &st_multi = *(Gudhi::Simplex_tree*)(splxptr); + int num_parameters = st_multi.get_number_of_parameters(); + std::vector out(degrees.size(), multi_filtration_grid(num_parameters)); + std::vector degree_index(degrees.size()); + int count = 0; + for (auto degree : degrees){ + degree_index[degree] = count; count++; + out[degree_index[degree]].reserve(st_multi.num_simplices()); + } + + for (const auto &simplex_handle : st_multi.complex_simplex_range()){ + const auto filtration = st_multi.filtration(simplex_handle); + const auto degree = st_multi.dimension(simplex_handle); + if (std::find(degrees.begin(), degrees.end(), degree) == degrees.end()) continue; + for (int parameter=0; parameter < num_parameters; parameter++){ + out[degree_index[degree]][parameter].push_back(filtration[parameter]); + } + } + return out; + +} + +} // namespace Gudhi + + + +namespace std { + +template<> +class numeric_limits +{ +public: + static Gudhi::multiparameter::multi_filtration_type infinity() throw(){ + return Gudhi::multiparameter::multi_filtration_type(1, std::numeric_limits::infinity()); + }; + + + static Gudhi::multiparameter::multi_filtration_type quiet_NaN() throw(){ + return Gudhi::multiparameter::multi_filtration_type(1, numeric_limits::quiet_NaN()); + }; + +}; + +} + + +#endif // SIMPLEX_TREE_MULTI_H_ 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 d7dbb541be..cb3205f816 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 @@ -41,7 +41,7 @@ struct GUDHI_EMPTY_BASE_CLASS_OPTIMIZATION Simplex_tree_node_explicit_storage : typedef typename SimplexTree::Simplex_key Simplex_key; Simplex_tree_node_explicit_storage(Siblings * sib = nullptr, - Filtration_value filtration = 0) + const Filtration_value& filtration = {}) : children_(sib) { this->assign_filtration(filtration); } diff --git a/src/Simplex_tree/include/gudhi/Simplex_tree/multi_filtrations/box.h b/src/Simplex_tree/include/gudhi/Simplex_tree/multi_filtrations/box.h new file mode 100644 index 0000000000..a4c28fdf1f --- /dev/null +++ b/src/Simplex_tree/include/gudhi/Simplex_tree/multi_filtrations/box.h @@ -0,0 +1,164 @@ +/* This file is part of the MMA Library - https://gitlab.inria.fr/dloiseau/multipers - which is released under MIT. + * See file LICENSE for full license details. + * Author(s): David Loiseaux + * + * Copyright (C) 2021 Inria + * + * Modification(s): + * + */ +/** + * @file box.h + * @author David Loiseaux + * @brief BOX. + */ + +#ifndef BOX_H_INCLUDED +#define BOX_H_INCLUDED + +#include +#include +#include +#include +#include + +#include "finitely_critical_filtrations.h" + + + + +/** + * @brief Simple box in $\mathbb R^n$ . + */ + +namespace Gudhi::multiparameter::multi_filtrations{ + +template +class Box +{ + + using point_type = Finitely_critical_multi_filtration; +public: + Box(); + Box(T a, T b, T c, T d) : bottomCorner_({a,b}), upperCorner_({c,d}) {}; + Box(const point_type& bottomCorner, const point_type& upperCorner); + Box(const std::pair& box); + + void inflate(T delta); + const point_type& get_bottom_corner() const; + const point_type& get_upper_corner() const; + point_type& get_bottom_corner(); + point_type& get_upper_corner(); + bool contains(const point_type& point) const; + void infer_from_filters(const std::vector &Filters_list); + bool is_trivial() const ; + std::pair get_pair() const{ + return {bottomCorner_,upperCorner_}; + } + std::pair get_pair(){ + return {bottomCorner_,upperCorner_}; + } + +private: + point_type bottomCorner_; + point_type upperCorner_; +}; + +template +inline Box::Box() +{} + +template +inline Box::Box(const point_type &bottomCorner, const point_type &upperCorner) + : bottomCorner_(bottomCorner), + upperCorner_(upperCorner) +{ + assert(bottomCorner.size() == upperCorner.size() + && bottomCorner <= upperCorner + && "This box is trivial !"); +} + +template +inline Box::Box(const std::pair &box) + : bottomCorner_(box.first), + upperCorner_(box.second) +{} + + +template +inline void Box::inflate(T delta) +{ + bottomCorner_ -= delta; + upperCorner_ += delta; +} + +template +inline void Box::infer_from_filters(const std::vector &Filters_list){ + int dimension = Filters_list.size(); + int nsplx = Filters_list[0].size(); + std::vector lower(dimension); + std::vector upper(dimension); + for (int i = 0; i < dimension; i++){ + T min = Filters_list[i][0]; + T max = Filters_list[i][0]; + for (int j=1; j +inline bool Box::is_trivial() const { + return bottomCorner_.empty() || upperCorner_.empty() || bottomCorner_.size() != upperCorner_.size(); +} + +template +inline const typename Box::point_type &Box::get_bottom_corner() const +{ + return bottomCorner_; +} + +template +inline const typename Box::point_type &Box::get_upper_corner() const +{ + return upperCorner_; +} + +template +inline typename Box::point_type &Box::get_bottom_corner() +{ + return bottomCorner_; +} + +template +inline typename Box::point_type &Box::get_upper_corner() +{ + return upperCorner_; +} + +template +inline bool Box::contains(const point_type &point) const +{ + if (point.size() != bottomCorner_.size()) return false; + + return bottomCorner_ <= point && point <= upperCorner_; +} + +template +std::ostream& operator<<(std::ostream& os, const Box& box) +{ + os << "Box -- Bottom corner : "; + os << box.get_bottom_corner(); + os << ", Top corner : "; + os << box.get_upper_corner(); + return os; +} + +} // namespace Gudhi + + +#endif // BOX_H_INCLUDED diff --git a/src/Simplex_tree/include/gudhi/Simplex_tree/multi_filtrations/finitely_critical_filtrations.h b/src/Simplex_tree/include/gudhi/Simplex_tree/multi_filtrations/finitely_critical_filtrations.h new file mode 100644 index 0000000000..a75c45c641 --- /dev/null +++ b/src/Simplex_tree/include/gudhi/Simplex_tree/multi_filtrations/finitely_critical_filtrations.h @@ -0,0 +1,150 @@ +#ifndef FINITELY_CRITICAL_FILTRATIONS_H_ +#define FINITELY_CRITICAL_FILTRATIONS_H_ + +#include +#include +#include + +namespace Gudhi::multiparameter::multi_filtrations{ + +template +class Finitely_critical_multi_filtration : public std::vector { + // Class to prevent doing illegal stuff with the standard library, e.g., compare two vectors +public: + Finitely_critical_multi_filtration() : std::vector() {}; + Finitely_critical_multi_filtration(int n) : std::vector(n, -std::numeric_limits::infinity()) {}; // minus infinity by default + Finitely_critical_multi_filtration(int n, T value) : std::vector(n,value) {}; + Finitely_critical_multi_filtration(std::initializer_list init) : std::vector(init) {}; + Finitely_critical_multi_filtration(const std::vector& v) : std::vector(v) {}; + Finitely_critical_multi_filtration(typename std::vector::iterator it_begin,typename std::vector::iterator it_end) : std::vector(it_begin, it_end) {}; + Finitely_critical_multi_filtration(typename std::vector::const_iterator it_begin,typename std::vector::const_iterator it_end) : std::vector(it_begin, it_end) {}; + + + operator std::vector&() const { + return *this; + } + std::vector get_vector() const{ + return static_cast>(*this); + } + + + + friend bool operator<(const Finitely_critical_multi_filtration& a, const Finitely_critical_multi_filtration& b) + { + bool isSame = true; + int n = std::min(a.size(), b.size()); + for (int i = 0; i < n; ++i){ + if (a[i] > b[i]) return false; + if (isSame && a[i] != b[i]) isSame = false; + } + if (isSame) return false; + return true; + } + friend bool operator<=(const Finitely_critical_multi_filtration& a, const Finitely_critical_multi_filtration& b) + { + int n = std::min(a.size(), b.size()); + for (int i = 0; i < n; ++i){ + if (a[i] > b[i]) return false; + } + return true; + } + + + + //GREATER THAN OPERATORS + friend bool operator>(const Finitely_critical_multi_filtration& a, const Finitely_critical_multi_filtration& b) + { + return b=(const Finitely_critical_multi_filtration& a, const Finitely_critical_multi_filtration& b) + { + return b<=a; + } + + Finitely_critical_multi_filtration& operator=(const Finitely_critical_multi_filtration& a){ + std::vector::operator=(a); + return *this; + } + + std::vector& _convert_back(){ + return *this; + } + + + + + friend Finitely_critical_multi_filtration& operator-=(Finitely_critical_multi_filtration &result, const Finitely_critical_multi_filtration &to_substract){ + std::transform(result.begin(), result.end(), to_substract.begin(),result.begin(), std::minus()); + return result; + } + friend Finitely_critical_multi_filtration& operator+=(Finitely_critical_multi_filtration &result, const Finitely_critical_multi_filtration &to_add){ + std::transform(result.begin(), result.end(), to_add.begin(),result.begin(), std::plus()); + return result; + } + + friend Finitely_critical_multi_filtration& operator-=(Finitely_critical_multi_filtration &result, const T &to_substract){ + // std::transform(result.begin(), result.end(), to_substract.begin(),result.begin(), std::minus()); + for (auto & truc : result){ + truc -= to_substract; + } + return result; + } + friend Finitely_critical_multi_filtration& operator+=(Finitely_critical_multi_filtration &result, const T &to_add){ + for (auto & truc : result){ + truc += to_add; + } + return result; + } + + static std::vector> to_python(const std::vector>& to_convert){ + return std::vector>(to_convert.begin(), to_convert.end()); + } + + + static std::vector> from_python(const std::vector>& to_convert){ + return std::vector>(to_convert.begin(), to_convert.end());; + } + void push_to(const Finitely_critical_multi_filtration& x){ + if (this->size() != x.size()){ + std::cerr << "Does only work with 1-critical filtrations ! Sizes " << this->size() << " and " << x.size() << "are different !" << std::endl; + throw std::logic_error("Bad sizes"); + } + for (unsigned int i = 0; i < x.size(); i++) + this->at(i) = this->at(i) > x[i] ? this->at(i) : x[i]; + } + // Warning, this function assumes that the comparisons checks have already been made ! + void insert_new(Finitely_critical_multi_filtration to_concatenate){ + this->insert( + this->end(), std::move_iterator(to_concatenate.begin()), std::move_iterator(to_concatenate.end()) + ); + } + + + // scalar product of a filtration value with x. + T linear_projection(const std::vector& x){ + T projection=0; + unsigned int size = std::min(x.size(), this->size()); + for (auto i =0u; iat(i); + return projection; + } + + // easy debug + friend std::ostream& operator<<(std::ostream& stream, const Finitely_critical_multi_filtration& truc){ + stream << "["; + for(unsigned int i = 0; i < truc.size()-1; i++){ + stream << truc[i] << ", "; + } + if(!truc.empty()) stream << truc.back(); + stream << "]"; + return stream; + } + + + + +}; + + +} // namespace Gudhi +#endif // FINITELY_CRITICAL_FILTRATIONS_H_ diff --git a/src/Simplex_tree/include/gudhi/Simplex_tree/multi_filtrations/line.h b/src/Simplex_tree/include/gudhi/Simplex_tree/multi_filtrations/line.h new file mode 100644 index 0000000000..0fa82cf715 --- /dev/null +++ b/src/Simplex_tree/include/gudhi/Simplex_tree/multi_filtrations/line.h @@ -0,0 +1,94 @@ +/* This file is part of the MMA Library - https://gitlab.inria.fr/dloiseau/multipers - which is released under MIT. + * See file LICENSE for full license details. + * Author(s): David Loiseaux + * + * Copyright (C) 2022 Inria + * + * Modification(s): + */ +/** + * @file line_filtration_translation.h + * @author David Loiseaux + * @brief + */ + + +#ifndef LINE_FILTRATION_TRANSLATION_H_INCLUDED +#define LINE_FILTRATION_TRANSLATION_H_INCLUDED + +#include "box.h" +#include "finitely_critical_filtrations.h" + +namespace Gudhi::multiparameter::multi_filtrations{ + + + template + class Line + { + + public: + using point_type = Finitely_critical_multi_filtration; + Line(); + Line(point_type x); + Line(point_type x, point_type v); + point_type push_forward(point_type x) const; + point_type push_back(point_type x) const; + int get_dim() const; + std::pair get_bounds(const Box &box) const; + + + private: + point_type basepoint_; // any point on the line + point_type direction_; // direction of the line + + }; + template + Line::Line(){} + + template + Line::Line(point_type x){ + this->basepoint_.swap(x); + // this->direction_ = {}; // diagonal line + } + template + Line::Line(point_type x, point_type v){ + this->basepoint_.swap(x); + this->direction_.swap(v); + } + template + typename Line::point_type Line::push_forward(point_type x) const { //TODO remove copy + x -= basepoint_; + T t = - std::numeric_limits::infinity();; + for (std::size_t i = 0; idirection_.size() > i ? direction_[i] : 1; + t = std::max(t, x[i]/dir); + } + point_type out(basepoint_.size()); + for (unsigned int i = 0; i < out.size(); i++) + out[i] = basepoint_[i] + t * (this->direction_.size() > i ? direction_[i] : 1) ; + return out; + } + template + typename Line::point_type Line::push_back(point_type x) const{ + x -= basepoint_; + T t = std::numeric_limits::infinity(); + for (unsigned int i = 0; idirection_.size() > i ? direction_[i] : 1; + t = std::min(t, x[i]/dir); + } + point_type out(basepoint_.size()); + for (unsigned int i = 0; i < out.size(); i++) + out[i] = basepoint_[i] + t * (this->direction_.size() > i ? direction_[i] : 1) ; + return out; + } + template + int Line::get_dim() const{ + return basepoint_.size(); + } + template + std::pair::point_type, typename Line::point_type> Line::get_bounds(const Box &box) const{ + return {this->push_forward(box.get_bottom_corner()), this->push_back(box.get_upper_corner())}; + } +} + +#endif // LINE_FILTRATION_TRANSLATION_H_INCLUDED 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 986f811b54..e119008670 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 @@ -36,6 +36,7 @@ struct Simplex_tree_options_stable_simplex_handles { 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_graph_expansion_unit_test.cpp b/src/Simplex_tree/test/simplex_tree_graph_expansion_unit_test.cpp index ed0f38d632..3a0eea485f 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_make_filtration_non_decreasing_unit_test.cpp b/src/Simplex_tree/test/simplex_tree_make_filtration_non_decreasing_unit_test.cpp index 05b7052612..3ac3b97401 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_unit_test.cpp b/src/Simplex_tree/test/simplex_tree_unit_test.cpp index fa1a444151..de1303e1e2 100644 --- a/src/Simplex_tree/test/simplex_tree_unit_test.cpp +++ b/src/Simplex_tree/test/simplex_tree_unit_test.cpp @@ -41,6 +41,7 @@ struct Simplex_tree_options_stable_simplex_handles { 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,