diff --git a/multipers/gudhi/Persistence_slices_interface.h b/multipers/gudhi/Persistence_slices_interface.h index 1ec86ab3..528c84a2 100644 --- a/multipers/gudhi/Persistence_slices_interface.h +++ b/multipers/gudhi/Persistence_slices_interface.h @@ -7,7 +7,8 @@ #include "gudhi/multi_simplex_tree_helpers.h" #include "gudhi/persistence_matrix_options.h" #include "gudhi/Multi_parameter_filtered_complex.h" -#include "gudhi/Multi_persistence/Persistence_interface_matrix.h" +#include "gudhi/Multi_persistence/Persistence_interface_vineyard.h" +#include "gudhi/Multi_persistence/Persistence_interface_homology.h" #include "gudhi/Multi_persistence/Persistence_interface_cohomology.h" #include "gudhi/Dynamic_multi_parameter_filtration.h" #include "gudhi/Degree_rips_bifiltration.h" @@ -39,26 +40,26 @@ enum Filtration_containers_strs : std::uint8_t { using Available_columns = Gudhi::persistence_matrix::Column_types; -template -struct Multi_persistence_options : Gudhi::persistence_matrix::Default_options { - using Index = std::uint32_t; - static const bool has_matrix_maximal_dimension_access = false; - static const bool has_column_pairings = true; - static const bool has_vine_update = true; - static const bool can_retrieve_representative_cycles = true; -}; - -template -struct Multi_persistence_Clement_options : Gudhi::persistence_matrix::Default_options { - using Index = std::uint32_t; - static const bool has_matrix_maximal_dimension_access = false; - static const bool has_column_pairings = true; - static const bool has_vine_update = true; - static const bool is_of_boundary_type = false; - static const Gudhi::persistence_matrix::Column_indexation_types column_indexation_type = - Gudhi::persistence_matrix::Column_indexation_types::POSITION; - static const bool can_retrieve_representative_cycles = true; -}; +// template +// struct Multi_persistence_options : Gudhi::persistence_matrix::Default_options { +// using Index = std::uint32_t; +// static const bool has_matrix_maximal_dimension_access = false; +// static const bool has_column_pairings = true; +// static const bool has_vine_update = true; +// static const bool can_retrieve_representative_cycles = true; +// }; + +// template +// struct Multi_persistence_Clement_options : Gudhi::persistence_matrix::Default_options { +// using Index = std::uint32_t; +// static const bool has_matrix_maximal_dimension_access = false; +// static const bool has_column_pairings = true; +// static const bool has_vine_update = true; +// static const bool is_of_boundary_type = false; +// static const Gudhi::persistence_matrix::Column_indexation_types column_indexation_type = +// Gudhi::persistence_matrix::Column_indexation_types::POSITION; +// static const bool can_retrieve_representative_cycles = true; +// }; template struct No_vine_multi_persistence_options : Gudhi::persistence_matrix::Default_options { @@ -68,50 +69,61 @@ struct No_vine_multi_persistence_options : Gudhi::persistence_matrix::Default_op static const bool has_vine_update = false; }; -template -struct fix_presentation_options : Gudhi::persistence_matrix::Default_options { - using Index = std::uint32_t; - static const bool has_row_access = row_access; - static const bool has_map_column_container = false; - static const bool has_removable_columns = false; // WARN : idx will change if map is not true +// template +// struct fix_presentation_options : Gudhi::persistence_matrix::Default_options { +// using Index = std::uint32_t; +// static const bool has_row_access = row_access; +// static const bool has_map_column_container = false; +// static const bool has_removable_columns = false; // WARN : idx will change if map is not true +// }; + +template +struct Multi_persistence_vineyard_ru_options : Gudhi::vineyard::Default_vineyard_options { + static constexpr bool is_RU = true; + static const Gudhi::persistence_matrix::Column_types column_type = col; }; -template -using BackendOptionsWithVine = Multi_persistence_options; +template +struct Multi_persistence_vineyard_chain_options : Gudhi::vineyard::Default_vineyard_options { + static constexpr bool is_RU = false; + static const Gudhi::persistence_matrix::Column_types column_type = col; +}; + +// template +// using BackendOptionsWithVine = Multi_persistence_options; template using BackendOptionsWithoutVine = No_vine_multi_persistence_options; -template -using ClementBackendOptionsWithVine = Multi_persistence_Clement_options; +// template +// using ClementBackendOptionsWithVine = Multi_persistence_Clement_options; // using SimplicialStructure = Gudhi::multiparameter::truc_interface::SimplicialStructure; template using StructureStuff = Gudhi::multi_persistence::Multi_parameter_filtered_complex; -template -using MatrixBackendNoVine = Gudhi::multi_persistence::Persistence_interface_matrix>; +template +using MatrixBackendNoVine = Gudhi::multi_persistence::Persistence_interface_homology, Filtration>; template -using MatrixBackendVine = Gudhi::multi_persistence::Persistence_interface_matrix>; +using MatrixBackendVine = Gudhi::multi_persistence::Persistence_interface_vineyard>; template -using ClementMatrixBackendVine = - Gudhi::multi_persistence::Persistence_interface_matrix>; +using ClementMatrixBackendVine = Gudhi::multi_persistence::Persistence_interface_vineyard>; template using GraphBackendVine = Gudhi::multiparameter::truc_interface::Persistence_backend_h0>; template using Filtration_value = Gudhi::multi_filtration::Multi_parameter_filtration; -template -using SimplicialNoVineMatrixTruc = Gudhi::multi_persistence::Slicer, MatrixBackendNoVine>; +template +using SimplicialNoVineMatrixTruc = Gudhi::multi_persistence::Slicer, MatrixBackendNoVine>; template using GeneralVineTruc = Gudhi::multi_persistence::Slicer, MatrixBackendVine>; -template -using GeneralNoVineTruc = Gudhi::multi_persistence::Slicer, MatrixBackendNoVine>; +template +using GeneralNoVineTruc = Gudhi::multi_persistence::Slicer, MatrixBackendNoVine>; template using GeneralVineClementTruc = Gudhi::multi_persistence::Slicer, ClementMatrixBackendVine>; @@ -128,8 +140,8 @@ using Multi_critical_filtration_value = Gudhi::multi_filtration::Multi_parameter template using KCriticalVineTruc = Gudhi::multi_persistence::Slicer, MatrixBackendVine>; -template -using Matrix_interface = std::conditional_t, MatrixBackendNoVine>; +template +using Matrix_interface = std::conditional_t, MatrixBackendNoVine>; template using filtration_options = std::conditional_t, Gudhi::multi_filtration::Degree_rips_bifiltration>>; -template -using MatrixTrucPythonInterface = - Gudhi::multi_persistence::Slicer, Matrix_interface>; +// template +// using MatrixTrucPythonInterface = Gudhi::multi_persistence::Slicer, Matrix_interface>; enum class BackendsEnum : std::uint8_t { Matrix, Graph, Clement, GudhiCohomology }; @@ -154,7 +165,7 @@ struct PersBackendOptsImpl; template struct PersBackendOptsImpl { - using type = Matrix_interface; + using type = Matrix_interface; }; template diff --git a/multipers/gudhi/gudhi/Degree_rips_bifiltration.h b/multipers/gudhi/gudhi/Degree_rips_bifiltration.h index 2fe87ab5..18a7cafe 100644 --- a/multipers/gudhi/gudhi/Degree_rips_bifiltration.h +++ b/multipers/gudhi/gudhi/Degree_rips_bifiltration.h @@ -31,13 +31,16 @@ #include #include +#ifdef GUDHI_USE_TBB +#include +#endif + #include #include #include #include #include #include -#include namespace Gudhi::multi_filtration { @@ -1708,7 +1711,8 @@ class Degree_rips_bifiltration GUDHI_CHECK_code(const OneDimArray &indices = grid[1]); const OneDimArray &values = grid[0]; - auto todo = [&](size_type g) { + + auto project_generator = [&](size_type g) { GUDHI_CHECK_code(GUDHI_CHECK(static_cast(indices[g]) == g, std::invalid_argument("Unvalid grid."))); auto v = static_cast(generators_[g]); @@ -1718,11 +1722,12 @@ class Degree_rips_bifiltration } generators_[g] = coordinate ? static_cast(d) : static_cast(values[d]); }; + #ifdef GUDHI_USE_TBB - tbb::parallel_for(size_type{0}, num_generators(), [&](size_type g) { todo(g); }); + tbb::parallel_for(size_type(0), num_generators(), project_generator); #else for (size_type g = 0; g < num_generators(); ++g) { - todo(g); + project_generator(g); } #endif } @@ -1842,9 +1847,7 @@ class Degree_rips_bifiltration #ifdef GUDHI_USE_TBB std::vector projections(f.num_generators()); - tbb::parallel_for(size_type{0}, f.num_generators(), [&](size_type g) { - projections[g] = project_generator(g); - }); + tbb::parallel_for(size_type{0}, f.num_generators(), [&](size_type g) { projections[g] = project_generator(g); }); if constexpr (Co) { return *std::max_element(projections.begin(), projections.end()); } else { diff --git a/multipers/gudhi/gudhi/Dynamic_multi_parameter_filtration.h b/multipers/gudhi/gudhi/Dynamic_multi_parameter_filtration.h index 24b4dd1e..602a54bb 100644 --- a/multipers/gudhi/gudhi/Dynamic_multi_parameter_filtration.h +++ b/multipers/gudhi/gudhi/Dynamic_multi_parameter_filtration.h @@ -31,10 +31,13 @@ #include #include +#ifdef GUDHI_USE_TBB +#include +#endif + #include #include #include -#include namespace Gudhi::multi_filtration { @@ -1856,9 +1859,8 @@ class Dynamic_multi_parameter_filtration std::invalid_argument("The grid should not be smaller than the number of parameters in the filtration value.")); #ifdef GUDHI_USE_TBB - tbb::parallel_for(size_type(0), num_generators(), [&](size_type i){ - generators_[i].project_onto_grid(grid, coordinate); - }); + tbb::parallel_for( + size_type(0), num_generators(), [&](size_type i) { generators_[i].project_onto_grid(grid, coordinate); }); #else for (Generator &g : generators_) { g.project_onto_grid(grid, coordinate); @@ -1946,6 +1948,17 @@ class Dynamic_multi_parameter_filtration { if (f.num_generators() == 1) return compute_linear_projection(f.generators_[0], x); +#ifdef GUDHI_USE_TBB + std::vector projections(f.num_generators()); + tbb::parallel_for(size_type{0}, f.num_generators(), [&](size_type g) { + projections[g] = compute_linear_projection(f.generators_[g], x); + }); + if constexpr (Co) { + return *std::max_element(projections.begin(), projections.end()); + } else { + return *std::min_element(projections.begin(), projections.end()); + } +#else if constexpr (Co) { U projection = std::numeric_limits::lowest(); for (const Generator &g : f.generators_) { @@ -1961,6 +1974,7 @@ class Dynamic_multi_parameter_filtration } return projection; } +#endif } /** diff --git a/multipers/gudhi/gudhi/Fields/Multi_field.h b/multipers/gudhi/gudhi/Fields/Multi_field.h index 3df57b9c..7fe0ac5b 100644 --- a/multipers/gudhi/gudhi/Fields/Multi_field.h +++ b/multipers/gudhi/gudhi/Fields/Multi_field.h @@ -2,7 +2,7 @@ * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. * Author(s): Hannah Schreiber, Clément Maria * - * Copyright (C) 2022-24 Inria + * Copyright (C) 2022 Inria * * Modification(s): * - YYYY/MM Author: Description of the modification @@ -305,7 +305,8 @@ class Multi_field_element const Characteristic& productOfCharacteristics) const { GUDHI_CHECK(productOfCharacteristics >= 0 && productOfCharacteristics <= productOfAllCharacteristics_, - "The given product is not the product of a subset of the current Multi-field characteristics."); + std::invalid_argument( + "The given product is not the product of a subset of the current Multi-field characteristics.")); Characteristic QR; mpz_gcd(QR.get_mpz_t(), element_.get_mpz_t(), productOfCharacteristics.get_mpz_t()); // QR <- gcd(x,QS) @@ -350,7 +351,8 @@ class Multi_field_element static Multi_field_element get_partial_multiplicative_identity(const Characteristic& productOfCharacteristics) { GUDHI_CHECK(productOfCharacteristics >= 0 && productOfCharacteristics <= productOfAllCharacteristics_, - "The given product is not the product of a subset of the current Multi-field characteristics."); + std::invalid_argument( + "The given product is not the product of a subset of the current Multi-field characteristics.")); if (productOfCharacteristics == 0 || productOfCharacteristics == productOfAllCharacteristics_) { return get_multiplicative_identity(); diff --git a/multipers/gudhi/gudhi/Fields/Multi_field_operators.h b/multipers/gudhi/gudhi/Fields/Multi_field_operators.h index 22dfd1cc..a5ec00d0 100644 --- a/multipers/gudhi/gudhi/Fields/Multi_field_operators.h +++ b/multipers/gudhi/gudhi/Fields/Multi_field_operators.h @@ -2,7 +2,7 @@ * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. * Author(s): Hannah Schreiber, Clément Maria * - * Copyright (C) 2022-24 Inria + * Copyright (C) 2022 Inria * * Modification(s): * - YYYY/MM Author: Description of the modification @@ -366,7 +366,8 @@ class Multi_field_operators const Characteristic& productOfCharacteristics) const { GUDHI_CHECK(productOfCharacteristics >= 0 && productOfCharacteristics <= productOfAllCharacteristics_, - "The given product is not the product of a subset of the current Multi-field characteristics."); + std::invalid_argument( + "The given product is not the product of a subset of the current Multi-field characteristics.")); Characteristic QR; mpz_gcd(QR.get_mpz_t(), e.get_mpz_t(), productOfCharacteristics.get_mpz_t()); // QR <- gcd(x,QS) @@ -409,7 +410,8 @@ class Multi_field_operators [[nodiscard]] Element get_partial_multiplicative_identity(const Characteristic& productOfCharacteristics) const { GUDHI_CHECK(productOfCharacteristics >= 0 && productOfCharacteristics <= productOfAllCharacteristics_, - "The given product is not the product of a subset of the current Multi-field characteristics."); + std::invalid_argument( + "The given product is not the product of a subset of the current Multi-field characteristics.")); if (productOfCharacteristics == nullCharacteristic || productOfCharacteristics == productOfAllCharacteristics_) { return get_multiplicative_identity(); diff --git a/multipers/gudhi/gudhi/Fields/Multi_field_shared.h b/multipers/gudhi/gudhi/Fields/Multi_field_shared.h index cb48c7b0..527229ad 100644 --- a/multipers/gudhi/gudhi/Fields/Multi_field_shared.h +++ b/multipers/gudhi/gudhi/Fields/Multi_field_shared.h @@ -2,7 +2,7 @@ * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. * Author(s): Hannah Schreiber, Clément Maria * - * Copyright (C) 2022-24 Inria + * Copyright (C) 2022 Inria * * Modification(s): * - YYYY/MM Author: Description of the modification @@ -347,7 +347,8 @@ class Shared_multi_field_element const Characteristic& productOfCharacteristics) const { GUDHI_CHECK(productOfCharacteristics >= 0 && productOfCharacteristics <= productOfAllCharacteristics_, - "The given product is not the product of a subset of the current Multi-field characteristics."); + std::invalid_argument( + "The given product is not the product of a subset of the current Multi-field characteristics.")); Element QR; mpz_gcd(QR.get_mpz_t(), element_.get_mpz_t(), productOfCharacteristics.get_mpz_t()); // QR <- gcd(x,QS) @@ -390,7 +391,8 @@ class Shared_multi_field_element static Shared_multi_field_element get_partial_multiplicative_identity(const Characteristic& productOfCharacteristics) { GUDHI_CHECK(productOfCharacteristics >= 0 && productOfCharacteristics <= productOfAllCharacteristics_, - "The given product is not the product of a subset of the current Multi-field characteristics."); + std::invalid_argument( + "The given product is not the product of a subset of the current Multi-field characteristics.")); if (productOfCharacteristics == 0 || productOfCharacteristics == productOfAllCharacteristics_) { return get_multiplicative_identity(); diff --git a/multipers/gudhi/gudhi/Fields/Multi_field_small.h b/multipers/gudhi/gudhi/Fields/Multi_field_small.h index 432eef1b..54ac4f52 100644 --- a/multipers/gudhi/gudhi/Fields/Multi_field_small.h +++ b/multipers/gudhi/gudhi/Fields/Multi_field_small.h @@ -2,7 +2,7 @@ * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. * Author(s): Hannah Schreiber, Clément Maria * - * Copyright (C) 2022-24 Inria + * Copyright (C) 2022 Inria * * Modification(s): * - YYYY/MM Author: Description of the modification @@ -353,7 +353,8 @@ class Multi_field_element_with_small_characteristics Characteristic productOfCharacteristics) const { GUDHI_CHECK(productOfCharacteristics >= 0 && productOfCharacteristics <= productOfAllCharacteristics_, - "The given product is not the product of a subset of the current Multi-field characteristics."); + std::invalid_argument( + "The given product is not the product of a subset of the current Multi-field characteristics.")); Characteristic gcd = std::gcd(element_, productOfAllCharacteristics_); @@ -401,7 +402,8 @@ class Multi_field_element_with_small_characteristics const Characteristic& productOfCharacteristics) { GUDHI_CHECK(productOfCharacteristics >= 0 && productOfCharacteristics <= productOfAllCharacteristics_, - "The given product is not the product of a subset of the current Multi-field characteristics."); + std::invalid_argument( + "The given product is not the product of a subset of the current Multi-field characteristics.")); if (productOfCharacteristics == 0 || productOfCharacteristics == productOfAllCharacteristics_) { return get_multiplicative_identity(); diff --git a/multipers/gudhi/gudhi/Fields/Multi_field_small_operators.h b/multipers/gudhi/gudhi/Fields/Multi_field_small_operators.h index 15e1534a..45785f43 100644 --- a/multipers/gudhi/gudhi/Fields/Multi_field_small_operators.h +++ b/multipers/gudhi/gudhi/Fields/Multi_field_small_operators.h @@ -2,7 +2,7 @@ * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. * Author(s): Hannah Schreiber, Clément Maria * - * Copyright (C) 2022-24 Inria + * Copyright (C) 2022 Inria * * Modification(s): * - YYYY/MM Author: Description of the modification @@ -327,7 +327,8 @@ class Multi_field_operators_with_small_characteristics const Characteristic& productOfCharacteristics) const { GUDHI_CHECK(productOfCharacteristics >= 0 && productOfCharacteristics <= productOfAllCharacteristics_, - "The given product is not the product of a subset of the current Multi-field characteristics."); + std::invalid_argument( + "The given product is not the product of a subset of the current Multi-field characteristics.")); Characteristic gcd = std::gcd(e, productOfAllCharacteristics_); @@ -368,7 +369,8 @@ class Multi_field_operators_with_small_characteristics [[nodiscard]] Element get_partial_multiplicative_identity(const Characteristic& productOfCharacteristics) const { GUDHI_CHECK(productOfCharacteristics >= 0 && productOfCharacteristics <= productOfAllCharacteristics_, - "The given product is not the product of a subset of the current Multi-field characteristics."); + std::invalid_argument( + "The given product is not the product of a subset of the current Multi-field characteristics.")); if (productOfCharacteristics == nullCharacteristic || productOfCharacteristics == productOfAllCharacteristics_) { return get_multiplicative_identity(); diff --git a/multipers/gudhi/gudhi/Fields/Multi_field_small_shared.h b/multipers/gudhi/gudhi/Fields/Multi_field_small_shared.h index a57a6d77..b46d2bfd 100644 --- a/multipers/gudhi/gudhi/Fields/Multi_field_small_shared.h +++ b/multipers/gudhi/gudhi/Fields/Multi_field_small_shared.h @@ -2,7 +2,7 @@ * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. * Author(s): Hannah Schreiber, Clément Maria * - * Copyright (C) 2022-24 Inria + * Copyright (C) 2022 Inria * * Modification(s): * - YYYY/MM Author: Description of the modification @@ -395,7 +395,8 @@ class Shared_multi_field_element_with_small_characteristics Characteristic productOfCharacteristics) const { GUDHI_CHECK(productOfCharacteristics >= 0 && productOfCharacteristics <= productOfAllCharacteristics_, - "The given product is not the product of a subset of the current Multi-field characteristics."); + std::invalid_argument( + "The given product is not the product of a subset of the current Multi-field characteristics.")); Characteristic gcd = std::gcd(element_, productOfAllCharacteristics_); @@ -443,7 +444,8 @@ class Shared_multi_field_element_with_small_characteristics const Characteristic& productOfCharacteristics) { GUDHI_CHECK(productOfCharacteristics >= 0 && productOfCharacteristics <= productOfAllCharacteristics_, - "The given product is not the product of a subset of the current Multi-field characteristics."); + std::invalid_argument( + "The given product is not the product of a subset of the current Multi-field characteristics.")); if (productOfCharacteristics == 0 || productOfCharacteristics == productOfAllCharacteristics_) { return Shared_multi_field_element_with_small_characteristics(multiplicativeID_); diff --git a/multipers/gudhi/gudhi/Fields/Z2_field.h b/multipers/gudhi/gudhi/Fields/Z2_field.h index 18e309d6..91a45170 100644 --- a/multipers/gudhi/gudhi/Fields/Z2_field.h +++ b/multipers/gudhi/gudhi/Fields/Z2_field.h @@ -2,7 +2,7 @@ * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. * Author(s): Hannah Schreiber * - * Copyright (C) 2022-24 Inria + * Copyright (C) 2022 Inria * * Modification(s): * - YYYY/MM Author: Description of the modification diff --git a/multipers/gudhi/gudhi/Fields/Zp_field.h b/multipers/gudhi/gudhi/Fields/Zp_field.h index 4998d4ec..a60b6665 100644 --- a/multipers/gudhi/gudhi/Fields/Zp_field.h +++ b/multipers/gudhi/gudhi/Fields/Zp_field.h @@ -2,7 +2,7 @@ * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. * Author(s): Hannah Schreiber * - * Copyright (C) 2022-24 Inria + * Copyright (C) 2022 Inria * * Modification(s): * - YYYY/MM Author: Description of the modification diff --git a/multipers/gudhi/gudhi/Matrix.h b/multipers/gudhi/gudhi/Matrix.h index 82d043e0..93b5baeb 100644 --- a/multipers/gudhi/gudhi/Matrix.h +++ b/multipers/gudhi/gudhi/Matrix.h @@ -2,7 +2,7 @@ * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. * Author(s): Hannah Schreiber * - * Copyright (C) 2022-24 Inria + * Copyright (C) 2022 Inria * * Modification(s): * - YYYY/MM Author: Description of the modification @@ -1412,13 +1412,13 @@ class Matrix /** * @brief Only available if @ref PersistenceMatrixOptions::can_retrieve_representative_cycles is true. Pre-computes * the representative cycles of the current state of the filtration represented by the matrix. It needs to be called - * before calling @ref get_representative_cycles if the matrix was modified since last call. Otherwise the old cycles - * will be returned. + * before calling @ref get_all_representative_cycles if the matrix was modified since last call. Otherwise the old + * cycles will be returned. * * @param dim If different from default value, only the cycles of the given dimension are updated. * All others are erased. */ - void update_representative_cycles(Dimension dim = get_null_value()); + void update_all_representative_cycles(Dimension dim = get_null_value()); /** * @brief Only available if @ref PersistenceMatrixOptions::can_retrieve_representative_cycles is true. Pre-computes * the representative cycle in the current matrix state of the given bar. It needs to be called @@ -1430,22 +1430,22 @@ class Matrix void update_representative_cycle(const Bar& bar); /** * @brief Only available if @ref PersistenceMatrixOptions::can_retrieve_representative_cycles is true. - * Returns all representative cycles of the current filtration. @ref update_representative_cycles has to be called + * Returns all representative cycles of the current filtration. @ref update_all_representative_cycles has to be called * first if a modification to the matrix has to be token into account since last call. * * @return A const reference to the vector of representative cycles. */ - const std::vector& get_representative_cycles(); + const std::vector& get_all_representative_cycles() const; /** * @brief Only available if @ref PersistenceMatrixOptions::can_retrieve_representative_cycles is true. - * Returns the cycle representing the given bar. @ref update_representative_cycles or + * Returns the cycle representing the given bar. @ref update_all_representative_cycles or * @ref update_representative_cycle have to be called first if a modification to the matrix has to be token into * account since last call. * * @param bar A bar from the current barcode. * @return A const reference to the cycle representing @p bar. */ - const Cycle& get_representative_cycle(const Bar& bar); + const Cycle& get_representative_cycle(const Bar& bar) const; private: using Underlying_matrix = std::conditional_t< @@ -2123,10 +2123,10 @@ inline typename Matrix::Index Matrix -inline void Matrix::update_representative_cycles(Dimension dim) +inline void Matrix::update_all_representative_cycles(Dimension dim) { static_assert(PersistenceMatrixOptions::can_retrieve_representative_cycles, "This method was not enabled."); - matrix_.update_representative_cycles(dim); + matrix_.update_all_representative_cycles(dim); } template @@ -2138,15 +2138,15 @@ inline void Matrix::update_representative_cycle(const template inline const std::vector::Cycle>& -Matrix::get_representative_cycles() +Matrix::get_all_representative_cycles() const { static_assert(PersistenceMatrixOptions::can_retrieve_representative_cycles, "This method was not enabled."); - return matrix_.get_representative_cycles(); + return matrix_.get_all_representative_cycles(); } template inline const typename Matrix::Cycle& -Matrix::get_representative_cycle(const Bar& bar) +Matrix::get_representative_cycle(const Bar& bar) const { static_assert(PersistenceMatrixOptions::can_retrieve_representative_cycles, "This method was not enabled."); return matrix_.get_representative_cycle(bar); diff --git a/multipers/gudhi/gudhi/Multi_filtration/Multi_parameter_generator.h b/multipers/gudhi/gudhi/Multi_filtration/Multi_parameter_generator.h index 9c6da19a..1eca0e1b 100644 --- a/multipers/gudhi/gudhi/Multi_filtration/Multi_parameter_generator.h +++ b/multipers/gudhi/gudhi/Multi_filtration/Multi_parameter_generator.h @@ -31,10 +31,13 @@ #include #include +#ifdef GUDHI_USE_TBB +#include +#endif + #include #include #include -#include namespace Gudhi::multi_filtration { @@ -1359,7 +1362,8 @@ class Multi_parameter_generator * 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) { + void project_onto_grid(const std::vector &grid, bool coordinate = true) + { GUDHI_CHECK( grid.size() >= generator_.size(), std::invalid_argument("The grid should not be smaller than the number of parameters in the filtration value.")); @@ -1368,7 +1372,7 @@ class Multi_parameter_generator if (!is_finite()) generator_.resize(grid.size(), generator_[0]); - auto todo = [&](size_type p) { + auto project_parameter = [&](size_type p) { const auto &filtration = grid[p]; auto v = static_cast(generator_[p]); auto d = std::distance(filtration.begin(), std::lower_bound(filtration.begin(), filtration.end(), v)); @@ -1377,11 +1381,12 @@ class Multi_parameter_generator } generator_[p] = coordinate ? static_cast(d) : static_cast(filtration[d]); }; + #ifdef GUDHI_USE_TBB - tbb::parallel_for(size_type(0), generator_.size(), todo); + tbb::parallel_for(size_type(0), generator_.size(), project_parameter); #else for (size_type p = 0; p < generator_.size(); ++p) { - todo(p); + project_parameter(p); } #endif } diff --git a/multipers/gudhi/gudhi/Multi_parameter_filtered_complex.h b/multipers/gudhi/gudhi/Multi_parameter_filtered_complex.h index a53201c7..d5af9541 100644 --- a/multipers/gudhi/gudhi/Multi_parameter_filtered_complex.h +++ b/multipers/gudhi/gudhi/Multi_parameter_filtered_complex.h @@ -26,10 +26,13 @@ #include #include +#ifdef GUDHI_USE_TBB +#include +#endif + #include #include //for lex order #include -#include namespace Gudhi { namespace multi_persistence { @@ -42,18 +45,20 @@ namespace multi_persistence { * @brief Class storing the boundaries, the dimensions and the filtration values of all cells composing a complex. * * @tparam MultiFiltrationValue Filtration value class respecting the @ref MultiFiltrationValue concept. + * @tparam I Index type. Default value: std::uint32_t. + * @tparam D Dimension type. Default value: int. */ -template +template class Multi_parameter_filtered_complex { public: - using Index = std::uint32_t; /**< Complex index type. */ + using Index = I; /**< Complex index type. */ using Filtration_value = MultiFiltrationValue; /**< Filtration value type. */ using T = typename Filtration_value::value_type; /**< Numerical type of an element in a filtration value. */ using Filtration_value_container = std::vector; /**< Filtration value container type. */ - using Boundary = std::vector; /**< Cell boundary type, represented by the complex indices of its faces. */ + using Boundary = std::vector; /**< Cell boundary type, represented by the complex indices of its faces. */ using Boundary_container = std::vector; /**< Boundary container type. */ - using Dimension = int; /**< Dimension type. */ + using Dimension = D; /**< Dimension type. */ using Dimension_container = std::vector; /**< Dimension container type. */ /** @@ -117,10 +122,10 @@ class Multi_parameter_filtered_complex /** * @brief Copy constructor. */ - template - Multi_parameter_filtered_complex(const Multi_parameter_filtered_complex& complex) - : boundaries_(complex.get_boundaries()), - dimensions_(complex.get_dimensions()), + template + Multi_parameter_filtered_complex(const Multi_parameter_filtered_complex& complex) + : boundaries_(complex.get_boundaries().begin(), complex.get_boundaries().end()), + dimensions_(complex.get_dimensions().begin(), complex.get_dimensions().end()), filtrationValues_(complex.get_filtration_values().size()), maxDimension_(complex.get_max_dimension()), isOrderedByDimension_(complex.is_ordered_by_dimension()) @@ -149,11 +154,12 @@ class Multi_parameter_filtered_complex /** * @brief Assign operator. */ - template - Multi_parameter_filtered_complex& operator=(const Multi_parameter_filtered_complex& other) + template + Multi_parameter_filtered_complex& operator=( + const Multi_parameter_filtered_complex& other) { - boundaries_ = other.get_boundaries(); - dimensions_ = other.get_dimensions(); + boundaries_ = Boundary_container(other.get_boundaries().begin(), other.get_boundaries().end()); + dimensions_ = Dimension_container(other.get_dimensions().begin(), other.get_dimensions().end()); const auto& fils = other.get_filtration_values(); filtrationValues_ = Filtration_value_container(fils.size()); for (Index i = 0; i < filtrationValues_.size(); ++i) { @@ -318,14 +324,17 @@ class Multi_parameter_filtered_complex */ void coarsen_on_grid(const std::vector >& grid, bool coordinate = true) { - // for (auto gen = 0U; gen < filtrationValues_.size(); ++gen) { - // filtrationValues_[gen].project_onto_grid(grid, coordinate); - // } +#ifdef GUDHI_USE_TBB tbb::parallel_for(Index(0), Index(filtrationValues_.size()), [&](Index gen) { // TODO : preallocate for tbb + // preallocate what? filtrationValues_[gen].project_onto_grid(grid, coordinate); }); - +#else + for (auto& fil : filtrationValues_) { + fil.project_onto_grid(grid, coordinate); + } +#endif } /** @@ -394,7 +403,7 @@ class Multi_parameter_filtered_complex { using namespace Gudhi::multi_filtration; using Return_filtration_value = decltype(std::declval().template as_type()); - using Return_complex = Multi_parameter_filtered_complex; + using Return_complex = Multi_parameter_filtered_complex; typename Return_complex::Filtration_value_container coords(complex.get_number_of_cycle_generators()); for (Index gen = 0U; gen < coords.size(); ++gen) { @@ -403,27 +412,6 @@ class Multi_parameter_filtered_complex return Return_complex(complex.boundaries_, complex.dimensions_, coords); } - // /** - // * @brief Compares two boundaries and returns true if and only if the size of the first is strictly smaller than - // the - // * second, or, if the two sizes are the same, the first is lexicographically strictly smaller than the second. - // * The boundaries are assumed to be ordered by increasing values. - // */ - // static bool boundary_is_strictly_smaller_than(const Boundary& b1, const Boundary& b2) { - // // we want faces to be smaller than proper cofaces - // if (b1.size() < b2.size()) return true; - // if (b1.size() > b2.size()) return false; - - // // lexico for others - // for (Index i = 0; i < b2.size(); ++i){ - // if (b1[i] < b2[i]) return true; - // if (b1[i] > b2[i]) return false; - // } - - // // equal - // return false; - // } - /** * @brief Outstream operator. */ diff --git a/multipers/gudhi/gudhi/Multi_parameter_filtration.h b/multipers/gudhi/gudhi/Multi_parameter_filtration.h index 0338e33c..8e2fd206 100644 --- a/multipers/gudhi/gudhi/Multi_parameter_filtration.h +++ b/multipers/gudhi/gudhi/Multi_parameter_filtration.h @@ -32,10 +32,13 @@ #include #include +#ifdef GUDHI_USE_TBB +#include +#endif + #include #include #include -#include namespace Gudhi::multi_filtration { @@ -1785,36 +1788,34 @@ class Multi_parameter_filtration val = coordinate ? static_cast(d) : static_cast(filtration[d]); }; +#ifdef GUDHI_USE_TBB // TODO: verify if this really makes a differences in the 1-critical case, otherwise just keep the general case if constexpr (Ensure1Criticality) { -#ifdef GUDHI_USE_TBB - tbb::parallel_for(size_type(0), num_parameters(), [&](size_type p){ - project_generator_value(generators_[p], grid[p]); - }); -#else - for (size_type p = 0; p < num_parameters(); ++p) { - project_generator_value(generators_[p], grid[p]); - } -#endif + tbb::parallel_for( + size_type(0), num_parameters(), [&](size_type p) { project_generator_value(generators_[p], grid[p]); }); } else { - -#ifdef GUDHI_USE_TBB - tbb::parallel_for(size_type(0), num_generators(), [&](size_type g) { tbb::parallel_for(size_type(0), num_parameters(), [&](size_type p) { project_generator_value(generator_view_(g, p), grid[p]); }); }); + if (!coordinate && num_generators() > 1) simplify(); + } #else + // TODO: verify if this really makes a differences in the 1-critical case, otherwise just keep the general case + if constexpr (Ensure1Criticality) { + for (size_type p = 0; p < num_parameters(); ++p) { + project_generator_value(generators_[p], grid[p]); + } + } else { for (size_type g = 0; g < num_generators(); ++g) { for (size_type p = 0; p < num_parameters(); ++p) { project_generator_value(generator_view_(g, p), grid[p]); } } -#endif - if (!coordinate && num_generators() > 1) simplify(); } +#endif } // FONCTIONNALITIES @@ -1892,6 +1893,15 @@ class Multi_parameter_filtration if (f.num_generators() == 1) return project_generator(0); +#ifdef GUDHI_USE_TBB + std::vector projections(f.num_generators()); + tbb::parallel_for(size_type{0}, f.num_generators(), [&](size_type g) { projections[g] = project_generator(g); }); + if constexpr (Co) { + return *std::max_element(projections.begin(), projections.end()); + } else { + return *std::min_element(projections.begin(), projections.end()); + } +#else if constexpr (Co) { U projection = std::numeric_limits::lowest(); for (size_type g = 0; g < f.num_generators(); ++g) { @@ -1907,6 +1917,7 @@ class Multi_parameter_filtration } return projection; } +#endif } /** diff --git a/multipers/gudhi/gudhi/Multi_persistence/Line.h b/multipers/gudhi/gudhi/Multi_persistence/Line.h index 182a0698..23834f5b 100644 --- a/multipers/gudhi/gudhi/Multi_persistence/Line.h +++ b/multipers/gudhi/gudhi/Multi_persistence/Line.h @@ -19,7 +19,6 @@ #ifndef MP_LINE_FILTRATION_H_INCLUDED #define MP_LINE_FILTRATION_H_INCLUDED -#include #include #include @@ -93,7 +92,7 @@ class Line Point_t x(basePoint_.size()); if (direction_.size() > 0) { - for (std::size_t i = 0; i < x.size(); i++) x[i] = basePoint_[i] + t * direction_[i]; + 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; @@ -275,6 +274,18 @@ class Line return {bottom, top}; } + /** + * @brief Outstream operator. + */ + friend std::ostream& operator<<(std::ostream& stream, const Line& line) + { + stream << "Line:\n"; + stream << "base point: " << line.basePoint_ << "\n"; + stream << "direction: " << line.direction_ << "\n"; + + return stream; + } + private: Point_t basePoint_; /**< Any point on the line. */ Point_t direction_; /**< Direction of the line. */ diff --git a/multipers/gudhi/gudhi/Multi_persistence/Multi_parameter_filtered_complex_pcoh_interface.h b/multipers/gudhi/gudhi/Multi_persistence/Multi_parameter_filtered_complex_pcoh_interface.h index 0cd1d09c..91784133 100644 --- a/multipers/gudhi/gudhi/Multi_persistence/Multi_parameter_filtered_complex_pcoh_interface.h +++ b/multipers/gudhi/gudhi/Multi_persistence/Multi_parameter_filtered_complex_pcoh_interface.h @@ -36,19 +36,21 @@ namespace multi_persistence { * @ref Multi_parameter_filtered_complex. * * @tparam MultiFiltrationValue Filtration value type used in @ref Multi_parameter_filtered_complex. + * @tparam I Index type. Default value: std::uint32_t. + * @tparam D Dimension type. Default value: int. */ -template +template class Multi_parameter_filtered_complex_pcoh_interface { public: - using Complex = Multi_parameter_filtered_complex; /**< Complex type */ - using Simplex_key = std::uint32_t; /**< Simplex_key type */ - using Simplex_handle = Simplex_key; /**< Simplex_handle type */ - using Filtration_value = Simplex_key; /**< Internal filtration value type */ - using Dimension = int; /**< Internal dimension type */ - using Map = std::vector; /**< Map type */ - using Filtration_simplex_range = Map; /**< Filtration_simplex_range type */ - using Boundary_simplex_range = Map; /**< Boundary_simplex_range type */ + using Complex = Multi_parameter_filtered_complex; /**< Complex type */ + using Simplex_key = typename Complex::Index; /**< Simplex_key type */ + using Simplex_handle = Simplex_key; /**< Simplex_handle type */ + using Filtration_value = Simplex_key; /**< Internal filtration value type */ + using Dimension = typename Complex::Dimension; /**< Internal dimension type */ + using Map = std::vector; /**< Map type */ + using Filtration_simplex_range = Map; /**< Filtration_simplex_range type */ + using Boundary_simplex_range = Map; /**< Boundary_simplex_range type */ /** * @brief Default constructor, storing null pointers.Can be tested with @ref is_initialized(), but any other method @@ -156,9 +158,9 @@ class Multi_parameter_filtered_complex_pcoh_interface ~Multi_parameter_filtered_complex_pcoh_interface() = default; Multi_parameter_filtered_complex_pcoh_interface &operator=( - const Multi_parameter_filtered_complex_pcoh_interface &other) = delete; + const Multi_parameter_filtered_complex_pcoh_interface &other) = default; Multi_parameter_filtered_complex_pcoh_interface &operator=( - Multi_parameter_filtered_complex_pcoh_interface &&other) noexcept = delete; + Multi_parameter_filtered_complex_pcoh_interface &&other) noexcept = default; /** * @brief Swap operator. @@ -171,24 +173,6 @@ class Multi_parameter_filtered_complex_pcoh_interface be1.keys_.swap(be2.keys_); } - /** - * @brief Reinitializes the interface with the new given complex and permutation. - * To use instead of the classical assign operator `operator=`. - */ - void reinitialize(const Complex &boundaries, const Map &permutation) - { - boundaries_ = &boundaries; - newToOldPerm_ = &permutation; - keys_ = Map(boundaries.get_number_of_cycle_generators(), -1); - } - - void reset() - { - boundaries_ = nullptr; - newToOldPerm_ = nullptr; - keys_.clear(); - } - /** * @brief Returns `true` if and only if all pointers are not null. */ @@ -208,6 +192,8 @@ class Multi_parameter_filtered_complex_pcoh_interface return sh == null_simplex() ? -1 : boundaries_->get_dimensions()[sh]; } + Complex const *get_complex_ptr() const { return boundaries_; } + // assumes that pcoh will assign the keys from 0 to n in order of filtrations void assign_key(Simplex_handle sh, Simplex_key key) { @@ -242,8 +228,7 @@ class Multi_parameter_filtered_complex_pcoh_interface return boundaries_->get_boundaries()[sh]; } - friend std::ostream &operator<<(std::ostream &stream, - const Multi_parameter_filtered_complex_pcoh_interface &complex) + friend std::ostream &operator<<(std::ostream &stream, const Multi_parameter_filtered_complex_pcoh_interface &complex) { stream << "[\n"; for (auto i : complex.filtration_simplex_range()) { diff --git a/multipers/gudhi/gudhi/Multi_persistence/Persistence_interface_cohomology.h b/multipers/gudhi/gudhi/Multi_persistence/Persistence_interface_cohomology.h index edf73a18..5dc1abb6 100644 --- a/multipers/gudhi/gudhi/Multi_persistence/Persistence_interface_cohomology.h +++ b/multipers/gudhi/gudhi/Multi_persistence/Persistence_interface_cohomology.h @@ -19,8 +19,6 @@ #include -#include - #include #include #include @@ -42,85 +40,113 @@ template class Persistence_interface_cohomology { public: + /** + * @brief Complex interface for @ref Gudhi::persistent_cohomology::Persistent_cohomology + */ using PCOH_complex = Multi_parameter_filtered_complex_pcoh_interface; using Complex = typename PCOH_complex::Complex; /**< Complex type */ using Dimension = typename PCOH_complex::Dimension; /**< Dimension type */ using Index = typename PCOH_complex::Simplex_key; /**< Index type */ - using Map = typename PCOH_complex::Map; /**< Map type */ using Bar = Gudhi::persistence_matrix::Persistence_interval; /**< Bar type */ - using Field_Zp = Gudhi::persistent_cohomology::Field_Zp; - using Persistent_cohomology = Gudhi::persistent_cohomology::Persistent_cohomology; + using Field_Zp = Gudhi::persistent_cohomology::Field_Zp; /**< Cohomology field. */ using Barcode = std::vector; /**< Barcode type */ - template + using Map = std::vector; /**< Permutation map for filtration order. */ + /** + * @brief Cohomology computation type + */ + using Persistent_cohomology = Gudhi::persistent_cohomology::Persistent_cohomology; + template using As_type = Persistence_interface_cohomology; /**< This type. */ - static constexpr const auto nullDeath = Bar::inf; - static constexpr const bool is_vine = false; /** False. */ - static constexpr const bool has_rep_cycles = false; /** False. */ + static constexpr const auto nullDeath = Bar::inf; /**< Death value of the barcode when the bar never died. */ + static constexpr const bool is_vine = false; /**< False: no cheap update. */ + static constexpr const bool has_rep_cycles = false; /**< False: no rep cycle methods enabled. */ - Persistence_interface_cohomology() : interface_(), barcode_() {} + Persistence_interface_cohomology() : order_(), interface_(), barcode_() {} - // `permutation` is assumed to have stable size, i.e., its address never changes - Persistence_interface_cohomology(const Complex& cpx, const Map& permutation) : interface_(cpx, permutation) + template + Persistence_interface_cohomology(const Complex& cpx, const Filtration_range& filtrationValues, bool ignoreInf = false) + : order_(cpx.get_number_of_cycle_generators()), interface_(cpx, order_), barcode_() { - _initialize(); + _initialize_order(cpx, filtrationValues, ignoreInf); + _initialize_persistence(); } - Persistence_interface_cohomology(const Persistence_interface_cohomology& other) = delete; - - // permutation is assumed to be the same than from the copied object, just its address can change - Persistence_interface_cohomology(const Persistence_interface_cohomology& other, const Map& permutation) - : interface_(other.interface_, permutation), barcode_(other.barcode_) + Persistence_interface_cohomology(const Persistence_interface_cohomology& other) + : order_(other.order_), interface_(other.interface_, order_), barcode_(other.barcode_) {} - Persistence_interface_cohomology(Persistence_interface_cohomology&& other) = delete; - - // permutation is assumed to be the same than from the moved object, just its address can change - Persistence_interface_cohomology(Persistence_interface_cohomology&& other, const Map& permutation) - : interface_(std::move(other.interface_), permutation), barcode_(std::move(other.barcode_)) + Persistence_interface_cohomology(Persistence_interface_cohomology&& other) noexcept + : order_(std::move(other.order_)), + interface_(std::move(other.interface_), order_), + barcode_(std::move(other.barcode_)) {} ~Persistence_interface_cohomology() = default; - Persistence_interface_cohomology& operator=(const Persistence_interface_cohomology& other) = delete; - Persistence_interface_cohomology& operator=(Persistence_interface_cohomology&& other) noexcept = delete; + Persistence_interface_cohomology& operator=(const Persistence_interface_cohomology& other) + { + order_ = other.order_; + interface_ = PCOH_complex(other.interface_, order_); + barcode_ = other.barcode_; + return *this; + } + + Persistence_interface_cohomology& operator=(Persistence_interface_cohomology&& other) noexcept + { + order_ = std::move(other.order_); + interface_ = PCOH_complex(std::move(other.interface_), order_); + barcode_ = std::move(other.barcode_); + return *this; + } // TODO: swap? - template - void reinitialize(const Complex& cpx, const Map& permutation) + /** + * @brief The `ignoreInf` argument is not ignored for this class. + */ + template + void initialize(const Complex& cpx, const Filtration_range& filtrationValues, bool ignoreInf = false) { - interface_.reinitialize(cpx, permutation); - _initialize(); + _initialize_order(cpx, filtrationValues, ignoreInf); // may trigger a relocation for order_ + interface_ = PCOH_complex(cpx, order_); + _initialize_persistence(); } - void reset() + /** + * @brief The `ignoreInf` argument is not ignored for this class. + */ + template + void update(const Filtration_range& filtrationValues, bool ignoreInf = false) { - interface_.reset(); - barcode_.clear(); + GUDHI_CHECK(is_initialized(), std::logic_error("Barcode can not be updated uninitialized.")); + + auto const* cpx = interface_.get_complex_ptr(); + + // should not trigger a relocation for order_ as complex size did not change + _initialize_order(*cpx, filtrationValues, ignoreInf); + _initialize_persistence(); } [[nodiscard]] bool is_initialized() const { return interface_.is_initialized(); } - Dimension get_dimension(Index i) const - { - GUDHI_CHECK(is_initialized(), "Dimension can not be computed uninitialized."); - return interface_.dimension(interface_.simplex(i)); - } + const Map& get_current_order() const { return order_; } const Barcode& get_barcode() { - GUDHI_CHECK(is_initialized(), "Barcode can not be computed uninitialized."); + GUDHI_CHECK(is_initialized(), std::logic_error("Barcode can not be computed uninitialized.")); return barcode_; } - /** - * @brief Outstream operator. - */ friend std::ostream& operator<<(std::ostream& stream, const Persistence_interface_cohomology& pers) { stream << "Complex:\n"; stream << pers.interface_ << "\n"; + stream << "Permutation:\n"; + for (auto v : pers.order_) { + stream << v << " "; + } + stream << "\n"; stream << "Barcode:\n"; for (const auto bar : pers.barcode_) { stream << bar << "\n"; @@ -131,10 +157,36 @@ class Persistence_interface_cohomology } private: + Map order_; PCOH_complex interface_; Barcode barcode_; - void _initialize() + template + void _initialize_order(const Complex& complex, const Filtration_range& filtrationValues, bool ignoreInf) + { + const auto& dimensions = complex.get_dimensions(); + order_.resize(complex.get_number_of_cycle_generators()); + std::iota(order_.begin(), order_.end(), 0); + std::sort(order_.begin(), order_.end(), [&](Index i, Index j) { + if (ignoreInf) { + if (filtrationValues[i] != MultiFiltrationValue::T_inf && filtrationValues[j] == MultiFiltrationValue::T_inf) + return true; + // all elements at inf are considered equal + if (filtrationValues[i] == MultiFiltrationValue::T_inf) return false; + } + if (dimensions[i] > dimensions[j]) return false; + if (dimensions[i] < dimensions[j]) return true; + // if filtration values are equal, we don't care about order, so considered the same object + return filtrationValues[i] < filtrationValues[j]; + }); + if (ignoreInf) { + Index end = order_.size(); + while (end > 0 && filtrationValues[order_[end - 1]] == MultiFiltrationValue::T_inf) --end; + order_.resize(end); + } + } + + void _initialize_persistence() { Persistent_cohomology pcoh(interface_, true); pcoh.init_coefficients(2); diff --git a/multipers/gudhi/gudhi/Multi_persistence/Persistence_interface_homology.h b/multipers/gudhi/gudhi/Multi_persistence/Persistence_interface_homology.h new file mode 100644 index 00000000..81f46f97 --- /dev/null +++ b/multipers/gudhi/gudhi/Multi_persistence/Persistence_interface_homology.h @@ -0,0 +1,255 @@ +/* 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) 2025 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +/** + * @file Persistence_interface_homology.h + * @author Hannah Schreiber + * @brief Contains the @ref Gudhi::multi_persistence::Persistence_interface_homology class. + */ + +#ifndef MP_PERSISTENCE_INTERFACE_HOMOLOGY_H_INCLUDED +#define MP_PERSISTENCE_INTERFACE_HOMOLOGY_H_INCLUDED + +#include +#include +#include + +#include + +#include +#include +#include + +namespace Gudhi { +namespace multi_persistence { + +/** + * @class Persistence_interface_homology Persistence_interface_homology.h \ + * gudhi/Multi_persistence/Persistence_interface_homology.h + * @ingroup multi_persistence + * + * @brief Interface respecting the @ref PersistenceAlgorithm concept to use @ref Slicer with the homology + * and representative cycle algorithms implemented in @ref Gudhi::persistence_matrix::Matrix. + * + * @tparam PosIdxPersistenceMatrixOptions Options respecting the + * @ref Gudhi::persistence_matrix::PersistenceMatrixOptions concept such that + * @ref Gudhi::persistence_matrix::PersistenceMatrixOptions::has_column_pairings is true. Is it also recommended + * that @ref Gudhi::persistence_matrix::PersistenceMatrixOptions::Index and + * @ref Gudhi::persistence_matrix::PersistenceMatrixOptions::Dimension correspond to the respective types in the used + * @ref Slicer. + * @tparam MultiFiltrationValue Value of @ref Slicer::Filtration_value. + */ +template +class Persistence_interface_homology +{ + public: + using Options = PosIdxPersistenceMatrixOptions; /**< Matrix options */ + using Dimension = typename Options::Dimension; /**< Dimension type */ + using Index = typename Options::Index; /**< Index type */ + using Complex = Multi_parameter_filtered_complex; /**< Complex type */ + using Matrix = Gudhi::persistence_matrix::Matrix; /**< Matrix type */ + using Map = std::vector; /**< Permutation map for filtration order */ + using Bar = typename Matrix::Bar; /**< Bar type */ + using Cycle = typename Matrix::Cycle; /**< Cycle type */ + template + using As_type = Persistence_interface_homology; /**< This type. */ + + static constexpr auto nullDeath = Bar::inf; /**< Death value of the barcode when the bar never died. */ + static constexpr auto nullDimension = Matrix::template get_null_value(); /**< None value for dimension. */ + static constexpr bool is_vine = false; /**< False: no cheap update. */ + /** + * @brief True if and only if PosIdxPersistenceMatrixOptions::can_retrieve_representative_cycles is true. + */ + static constexpr bool has_rep_cycles = Options::can_retrieve_representative_cycles; + + Persistence_interface_homology() : complex_(nullptr) {} + + template + Persistence_interface_homology(const Complex& cpx, const Filtration_range& filtrationValues, bool ignoreInf = false) + : complex_(&cpx), matrix_(filtrationValues.size()), order_(filtrationValues.size()) + { + static_assert(Options::has_column_pairings, "Underlying matrix has to store barcode."); + + _initialize_order(filtrationValues, ignoreInf); + _initialize_persistence(); + } + + Persistence_interface_homology(const Persistence_interface_homology& other) = default; + + Persistence_interface_homology(Persistence_interface_homology&& other) noexcept + : complex_(std::exchange(other.complex_, nullptr)), + matrix_(std::move(other.matrix_)), + order_(std::move(other.order_)) + {} + + ~Persistence_interface_homology() = default; + + Persistence_interface_homology& operator=(const Persistence_interface_homology& other) = default; + + Persistence_interface_homology& operator=(Persistence_interface_homology&& other) noexcept + { + complex_ = std::exchange(other.complex_, nullptr); + matrix_ = std::move(other.matrix_); + order_ = std::move(other.order_); + return *this; + } + + // TODO: swap? + + /** + * @brief The `ignoreInf` argument is not ignored for this class. + */ + template + void initialize(const Complex& cpx, const Filtration_range& filtrationValues, bool ignoreInf = false) + { + complex_ = &cpx; + matrix_ = Matrix(filtrationValues.size()); + order_.resize(filtrationValues.size()); + + _initialize_order(filtrationValues, ignoreInf); + _initialize_persistence(); + } + + /** + * @brief The `ignoreInf` argument is not ignored for this class. + */ + template + void update(const Filtration_range& filtrationValues, bool ignoreInf = false) + { + GUDHI_CHECK(is_initialized(), std::logic_error("Barcode can not be updated uninitialized.")); + + initialize(*complex_, filtrationValues, ignoreInf); + } + + [[nodiscard]] bool is_initialized() const { return complex_ != nullptr; } + + const Map& get_current_order() const { return order_; } + + auto get_barcode() + { + GUDHI_CHECK(is_initialized(), std::logic_error("Barcode can not be computed uninitialized.")); + + const auto& barcode = matrix_.get_current_barcode(); + return boost::adaptors::transform(barcode, [&](const Bar& bar) -> Bar { + return Bar(order_[bar.birth], bar.death == Bar::inf ? bar.death : order_[bar.death], bar.dim); + }); + } + + auto get_all_representative_cycles(bool update = true, Dimension dim = nullDimension) + { + static_assert(has_rep_cycles, "`get_all_representative_cycles` is not enabled with the given options."); + GUDHI_CHECK(is_initialized(), std::logic_error("Representative cycles can not be computed uninitialized.")); + + if (update) matrix_.update_all_representative_cycles(dim); + const auto& cycles = matrix_.get_all_representative_cycles(); + return boost::adaptors::transform(cycles, [&](const Cycle& cycle) -> Cycle { + Cycle c(cycle.size()); + for (Index i = 0; i < cycle.size(); ++i) { + c[i] = order_[cycle[i]]; + } + return c; + }); + } + + auto get_representative_cycle(Index barcodeIndex, bool update = true) + { + static_assert(has_rep_cycles, "`get_representative_cycle` is not enabled with the given options."); + GUDHI_CHECK(is_initialized(), std::logic_error("Representative cycles can not be computed uninitialized.")); + + const auto& bar = matrix_.get_current_barcode()[barcodeIndex]; + if (update) matrix_.update_representative_cycle(bar); + const auto& cycle = matrix_.get_representative_cycle(bar); + return boost::adaptors::transform(cycle, [&](const Index& i) -> Index { return order_[i]; }); + } + + friend std::ostream& operator<<(std::ostream& stream, Persistence_interface_homology& pers) + { + stream << "Matrix:\n"; + stream << "[\n"; + for (auto i = 0U; i < pers.matrix_.get_number_of_columns(); i++) { + stream << "["; + for (const auto& v : pers.matrix_.get_column(i)) stream << v << ", "; + stream << "]\n"; + } + stream << "]\n"; + stream << "Permutation:\n"; + for (auto v : *pers.order_) { + stream << v << " "; + } + stream << "\n"; + + return stream; + } + + private: + Complex const* complex_; /**< Pointer to complex. */ + Matrix matrix_; + Map order_; + + template + void _initialize_order(const Filtration_range& filtrationValues, bool ignoreInf) + { + // GUDHI_CHECK_code(const auto& boundaryMatrix = complex_->get_boundaries();) + // Cannot use GUDHI_CHECK_code here because of the added GUDHI_ASSUME in multipers + // TODO: Things need to be reorganize later + [[maybe_unused]] const auto& boundaryMatrix = complex_->get_boundaries(); + const auto& dimensions = complex_->get_dimensions(); + + GUDHI_CHECK(boundaryMatrix.size() == dimensions.size(), + std::invalid_argument("Boundary and dimension range sizes are not matching.")); + GUDHI_CHECK(boundaryMatrix.size() == filtrationValues.size(), + std::invalid_argument("Boundary and filtration value range sizes are not matching.")); + + std::iota(order_.begin(), order_.end(), 0); + std::sort(order_.begin(), order_.end(), [&](Index i, Index j) { + if (ignoreInf) { + if (filtrationValues[i] != MultiFiltrationValue::T_inf && filtrationValues[j] == MultiFiltrationValue::T_inf) + return true; + // all elements at inf are considered equal + if (filtrationValues[i] == MultiFiltrationValue::T_inf) return false; + } + if (dimensions[i] < dimensions[j]) return true; + if (dimensions[i] > dimensions[j]) return false; + return filtrationValues[i] < filtrationValues[j]; + }); + if (ignoreInf) { + Index end = order_.size(); + while (end > 0 && filtrationValues[order_[end - 1]] == MultiFiltrationValue::T_inf) --end; + order_.resize(end); + } + } + + void _initialize_persistence() + { + const auto& boundaryMatrix = complex_->get_boundaries(); + const auto& dimensions = complex_->get_dimensions(); + + // simplex IDs need to be increasing in order, so the original ones cannot be used + Map orderInv(boundaryMatrix.size()); + Map translatedBoundary; + Index id = 0; + for (auto i : order_) { + orderInv[i] = id; // order is assumed to be a valid filtration + translatedBoundary.resize(boundaryMatrix[i].size()); + for (Index j = 0; j < boundaryMatrix[i].size(); ++j) { + translatedBoundary[j] = orderInv[boundaryMatrix[i][j]]; + } + std::sort(translatedBoundary.begin(), translatedBoundary.end()); + matrix_.insert_boundary(id, translatedBoundary, dimensions[i]); + ++id; // IDs corresponds to the indices in order_ + } + } +}; + +} // namespace multi_persistence +} // namespace Gudhi + +#endif // MP_PERSISTENCE_INTERFACE_HOMOLOGY_H_INCLUDED diff --git a/multipers/gudhi/gudhi/Multi_persistence/Persistence_interface_matrix.h b/multipers/gudhi/gudhi/Multi_persistence/Persistence_interface_matrix.h deleted file mode 100644 index be772506..00000000 --- a/multipers/gudhi/gudhi/Multi_persistence/Persistence_interface_matrix.h +++ /dev/null @@ -1,463 +0,0 @@ -/* 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) 2025 Inria - * - * Modification(s): - * - YYYY/MM Author: Description of the modification - */ - -/** - * @file Persistence_interface_matrix.h - * @author Hannah Schreiber - * @brief Contains the @ref Gudhi::multi_persistence::Persistence_interface_matrix class. - */ - -#ifndef MP_PERSISTENCE_INTERFACE_MATRIX_H_INCLUDED -#define MP_PERSISTENCE_INTERFACE_MATRIX_H_INCLUDED - -#include -#include -#include -#include - -#include - -#include -#include - -namespace Gudhi { -namespace multi_persistence { - -/** - * @class Persistence_interface_matrix Persistence_interface_matrix.h \ - * gudhi/Multi_persistence/Persistence_interface_matrix.h - * @ingroup multi_persistence - * - * @brief Interface respecting the @ref PersistenceAlgorithm concept to use @ref Slicer with the homology, vineyard - * and representative cycle algorithms implemented in @ref Gudhi::persistence_matrix::Matrix. - * - * @tparam PosIdxPersistenceMatrixOptions Options respecting the - * @ref Gudhi::persistence_matrix::PersistenceMatrixOptions concept such that, either - * @ref Gudhi::persistence_matrix::PersistenceMatrixOptions::column_indexation_type "column_indexation_type" is - * @ref Gudhi::persistence_matrix::Column_indexation_types::POSITION "POSITION", or, - * @ref Gudhi::persistence_matrix::PersistenceMatrixOptions::is_of_boundary_type "is_of_boundary_type" is true and - * @ref Gudhi::persistence_matrix::PersistenceMatrixOptions::column_indexation_type "column_indexation_type" is - * @ref Gudhi::persistence_matrix::Column_indexation_types::CONTAINER "CONTAINER". - */ -template -class Persistence_interface_matrix -{ - public: - using Options = PosIdxPersistenceMatrixOptions; - using Matrix = Gudhi::persistence_matrix::Matrix; /**< Complex type */ - using Dimension = typename Options::Dimension; /**< Dimension type */ - using Index = typename Options::Index; /**< Index type */ - using Map = std::vector; /**< Map type */ - using Bar = typename Matrix::Bar; /**< Bar type */ - using Cycle = typename Matrix::Cycle; /**< Cycle type */ - template - using As_type = Persistence_interface_matrix; /**< This type. */ - - class Barcode_iterator - : public boost::iterator_facade - { - private: - using Base = boost::iterator_facade; - - public: - using Barcode = typename Matrix::Barcode; - - using difference_type = typename Base::difference_type; - using size_type = std::size_t; - using const_reference = const Bar&; - - Barcode_iterator(const Barcode& barcode, const Map& permutation, size_type pos) - : barcode_(&barcode), - permutation_(&permutation), - currPos_(pos > barcode.size() ? barcode.size() : pos), - currBar_() - { - // end otherwise - if (currPos_ < barcode.size()) _update_bar(); - } - - // necessary for boost::iterator_range to be able to use operator(). - // operator[] not possible in this particular case - Bar operator[](difference_type n) const - { - const auto& b = (*barcode_)[currPos_ + n]; // n is out of range if currPos_ + n is out of range from barcode_ - return Bar(_posToIndex(b.birth), b.death == Bar::inf ? b.death : _posToIndex(b.death), b.dim); - } - - private: - using Pos_index = typename std::tuple_element<2, Bar>::type; - - friend class boost::iterator_core_access; - - bool equal(Barcode_iterator const& other) const - { - return barcode_ == other.barcode_ && permutation_ == other.permutation_ && currPos_ == other.currPos_; - } - - const_reference dereference() const { return currBar_; } - - void increment() - { - ++currPos_; - if (currPos_ < barcode_->size()) _update_bar(); - } - - void decrement() - { - --currPos_; - if (currPos_ < barcode_->size()) _update_bar(); - } - - void advance(difference_type n) - { - currPos_ += n; - if (currPos_ < barcode_->size()) _update_bar(); - } - - difference_type distance_to(const Barcode_iterator& other) const { return other.currPos_ - currPos_; } - - Pos_index _posToIndex(Pos_index pos) const { return (*permutation_)[pos]; } - - void _update_bar() - { - const auto& b = (*barcode_)[currPos_]; - currBar_.dim = b.dim; - currBar_.birth = _posToIndex(b.birth); - currBar_.death = b.death == Bar::inf ? b.death : _posToIndex(b.death); - } - - Barcode const* barcode_; - Map const* permutation_; - size_type currPos_; - Bar currBar_; - }; - - class Cycles_iterator - : public boost::iterator_facade - { - private: - using Base = boost::iterator_facade; - - public: - using Cycles = std::vector; - - using difference_type = typename Base::difference_type; - using size_type = std::size_t; - using const_reference = const Cycle&; - - Cycles_iterator(const Cycles& cycles, const Map& permutation, Map const* idToPos, size_type pos) - : cycles_(&cycles), - permutation_(&permutation), - idToPos_(idToPos), - currPos_(pos > cycles.size() ? cycles.size() : pos), - currCycle_() - { - if constexpr (Options::has_vine_update && !Options::is_of_boundary_type) { - GUDHI_CHECK(idToPos_ != nullptr, "ID to position map has to be set for chain matrices using vineyard."); - } - // end otherwise - if (currPos_ < cycles.size()) _update_cycle(); - } - - // necessary for boost::iterator_range to be able to use operator(). - // operator[] not possible in this particular case - Cycle operator[](difference_type n) const - { - Cycle res((*cycles_)[currPos_ + n]); - for (auto& id : res) { - if constexpr (PosIdxPersistenceMatrixOptions::is_z2) { - id = _id_to_index(id); - } else { - id.first = _id_to_index(id.first); - } - } - return res; - } - - private: - friend class boost::iterator_core_access; - - bool equal(Cycles_iterator const& other) const - { - return cycles_ == other.cycles_ && permutation_ == other.permutation_ && currPos_ == other.currPos_; - } - - const_reference dereference() const { return currCycle_; } - - void increment() - { - ++currPos_; - if (currPos_ < cycles_->size()) _update_cycle(); - } - - void decrement() - { - --currPos_; - if (currPos_ < cycles_->size()) _update_cycle(); - } - - void advance(difference_type n) - { - currPos_ += n; - if (currPos_ < cycles_->size()) _update_cycle(); - } - - difference_type distance_to(const Cycles_iterator& other) const { return other.currPos_ - currPos_; } - - Index _id_to_index(Index id) const - { - if constexpr (Options::is_of_boundary_type || !Options::has_vine_update) { - // works for RU because id == pos, but does not work for chain with vine - // we need a id to pos map in that case - return (*permutation_)[id]; - } else { - return (*permutation_)[(*idToPos_)[id]]; - } - } - - void _update_cycle() - { - const auto& c = (*cycles_)[currPos_]; - currCycle_.resize(c.size()); - for (size_type i = 0; i < c.size(); ++i) { - if constexpr (PosIdxPersistenceMatrixOptions::is_z2) { - currCycle_[i] = _id_to_index(c[i]); - } else { - currCycle_[i].first = _id_to_index(c[i]); - } - } - } - - Cycles const* cycles_; - Map const* permutation_; - Map const* idToPos_; - size_type currPos_; - Cycle currCycle_; - }; - - using Barcode = boost::iterator_range; /**< Barcode type */ - using Cycles = boost::iterator_range; /**< Cycle container type */ - - static constexpr const auto nullDeath = Bar::inf; - /** - * @brief True if and only if PosIdxPersistenceMatrixOptions::has_vine_update is true. - */ - static constexpr const bool is_vine = Options::has_vine_update; - /** - * @brief True if and only if PosIdxPersistenceMatrixOptions::can_retrieve_representative_cycles is true. - */ - static constexpr const bool has_rep_cycles = Options::can_retrieve_representative_cycles; - - Persistence_interface_matrix() : permutation_(nullptr) {} - - // `permutation` is assumed to have stable size, i.e., its address never changes - template - Persistence_interface_matrix(const Complex& cpx, const Map& permutation) - : matrix_(permutation.size()), permutation_(&permutation) - { - static_assert( - Options::column_indexation_type == Gudhi::persistence_matrix::Column_indexation_types::POSITION || - (Options::is_of_boundary_type && - Options::column_indexation_type == Gudhi::persistence_matrix::Column_indexation_types::CONTAINER), - "Matrix has a non supported index scheme."); - - _initialize(cpx); - } - - Persistence_interface_matrix(const Persistence_interface_matrix& other) = delete; - - // permutation is assumed to be the same than from the copied object, just its address can change - Persistence_interface_matrix(const Persistence_interface_matrix& other, const Map& permutation) - : matrix_(other.matrix_), permutation_(other.is_initialized() ? &permutation : nullptr), idToPos_(other.idToPos_) - { - GUDHI_CHECK(!other.is_initialized() || permutation == *other.permutation_, - "Only the address of the permutation vector is allowed to change, not its content."); - } - - Persistence_interface_matrix(Persistence_interface_matrix&& other) = delete; - - // permutation is assumed to be the same than from the moved object, just its address can change - Persistence_interface_matrix(Persistence_interface_matrix&& other, const Map& permutation) - : matrix_(std::move(other.matrix_)), - permutation_(other.is_initialized() ? &permutation : nullptr), - idToPos_(std::move(other.idToPos_)) - { - other.permutation_ = nullptr; - } - - ~Persistence_interface_matrix() = default; - - Persistence_interface_matrix& operator=(const Persistence_interface_matrix& other) = delete; - Persistence_interface_matrix& operator=(Persistence_interface_matrix&& other) noexcept = delete; - - // TODO: swap? - - template - void reinitialize(const Complex& cpx, const Map& permutation) - { - matrix_ = Matrix(permutation.size()); - permutation_ = &permutation; - _initialize(cpx); - } - - void reset() - { - matrix_ = Matrix(); - permutation_ = nullptr; - if constexpr (Options::has_vine_update && !Options::is_of_boundary_type) { - if (idToPos_) idToPos_->clear(); - } - } - - [[nodiscard]] bool is_initialized() const { return permutation_ != nullptr; } - - Dimension get_dimension(Index i) const - { - GUDHI_CHECK(is_initialized(), "Dimension can not be computed uninitialized."); - return matrix_.get_column_dimension(i); - } - - Barcode get_barcode() - { - GUDHI_CHECK(is_initialized(), "Barcode can not be computed uninitialized."); - - const auto& barcode = matrix_.get_current_barcode(); - return Barcode(Barcode_iterator(barcode, *permutation_, 0), - Barcode_iterator(barcode, *permutation_, barcode.size())); - } - - void vine_swap(Index i) - { - static_assert(is_vine, "`vine_swap` is not enabled with the given options."); - GUDHI_CHECK(is_initialized(), "Vineyard can not be computed uninitialized."); - - if constexpr (!Options::is_of_boundary_type) { - auto id1 = matrix_.get_pivot(i); - auto id2 = matrix_.get_pivot(i + 1); - std::swap((*idToPos_)[id1], (*idToPos_)[id2]); - } - matrix_.vine_swap(i); - } - - Cycles get_representative_cycles(bool update) - { - static_assert(has_rep_cycles, "`get_representative_cycles` is not enabled with the given options."); - GUDHI_CHECK(is_initialized(), "Representative cycles can not be computed uninitialized."); - - if (update) matrix_.update_representative_cycles(); - - const auto& cycles = matrix_.get_representative_cycles(); - Map* idToPosPtr = nullptr; - if constexpr (Options::has_vine_update && !Options::is_of_boundary_type) { - idToPosPtr = &(*idToPos_); - } - return Cycles(Cycles_iterator(cycles, *permutation_, idToPosPtr, 0), - Cycles_iterator(cycles, *permutation_, idToPosPtr, cycles.size())); - } - - Cycle get_representative_cycle(Index barcodeIndex, bool update) - { - static_assert(has_rep_cycles, "`get_representative_cycle` is not enabled with the given options."); - GUDHI_CHECK(is_initialized(), "Representative cycles can not be computed uninitialized."); - - auto id_to_index = [&](Index id) -> Index { - if constexpr (Options::is_of_boundary_type || !Options::has_vine_update) { - // works for RU because id == pos, but does not work for chain with vine - // we need a id to pos map in that case - return (*permutation_)[id]; - } else { - return (*permutation_)[(*idToPos_)[id]]; - } - }; - - if (update) matrix_.update_representative_cycles(); - - const auto& c = matrix_.get_representative_cycle(matrix_.get_current_barcode()[barcodeIndex]); - - Cycle cycle(c.size()); - for (Index i = 0; i < c.size(); ++i) { - if constexpr (PosIdxPersistenceMatrixOptions::is_z2) { - cycle[i] = id_to_index(c[i]); - } else { - cycle[i].first = id_to_index(c[i]); - } - } - return cycle; - } - - /** - * @brief Outstream operator. - */ - friend std::ostream& operator<<(std::ostream& stream, Persistence_interface_matrix& pers) - { - stream << "Matrix:\n"; - stream << "[\n"; - for (auto i = 0U; i < pers.matrix_.get_number_of_columns(); i++) { - stream << "["; - for (const auto& v : pers.matrix_.get_column(i)) stream << v << ", "; - stream << "]\n"; - } - stream << "]\n"; - stream << "Permutation:\n"; - if (pers.permutation_ != nullptr) { - for (auto v : *pers.permutation_) { - stream << v << " "; - } - } - stream << "\n"; - if constexpr (Options::has_vine_update && !Options::is_of_boundary_type) { - stream << "ID to position map:\n"; - for (auto v : *pers.idToPos_) { - stream << v << " "; - } - stream << "\n"; - } - - return stream; - } - - private: - Matrix matrix_; - Map const* permutation_; - std::optional idToPos_; - - template - void _initialize(const Complex& cpx) - { - if constexpr (Options::has_vine_update && !Options::is_of_boundary_type) { - idToPos_.emplace(); - idToPos_->reserve(permutation_->size()); - } - const auto& boundaries = cpx.get_boundaries(); - const auto& dimensions = cpx.get_dimensions(); - // simplex IDs need to be increasing in order, so the original ones cannot be used - Map permutationInv(cpx.get_number_of_cycle_generators()); - Map translated_boundary; - std::size_t id = 0; - for (auto i : *permutation_) { - permutationInv[i] = id; // permutation is assumed to be a valid filtration - translated_boundary.resize(boundaries[i].size()); - for (std::size_t j = 0; j < boundaries[i].size(); ++j) { - translated_boundary[j] = permutationInv[boundaries[i][j]]; - } - std::sort(translated_boundary.begin(), translated_boundary.end()); - matrix_.insert_boundary(id, translated_boundary, dimensions[i]); - if constexpr (Options::has_vine_update && !Options::is_of_boundary_type) { - idToPos_->push_back(id); - } - ++id; // IDs corresponds to the indices in permutation_ - } - } -}; - -} // namespace multi_persistence -} // namespace Gudhi - -#endif // MP_PERSISTENCE_INTERFACE_MATRIX_H_INCLUDED diff --git a/multipers/gudhi/gudhi/Multi_persistence/Persistence_interface_vineyard.h b/multipers/gudhi/gudhi/Multi_persistence/Persistence_interface_vineyard.h new file mode 100644 index 00000000..b6644ff5 --- /dev/null +++ b/multipers/gudhi/gudhi/Multi_persistence/Persistence_interface_vineyard.h @@ -0,0 +1,138 @@ +/* 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) 2026 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +/** + * @file Persistence_interface_vineyard.h + * @author Hannah Schreiber + * @brief Contains the @ref Gudhi::multi_persistence::Persistence_interface_vineyard class. + */ + +#ifndef MP_PERSISTENCE_INTERFACE_VINEYARD_H_INCLUDED +#define MP_PERSISTENCE_INTERFACE_VINEYARD_H_INCLUDED + +#include + +#include +#include + +namespace Gudhi { +namespace multi_persistence { + +/** + * @class Persistence_interface_vineyard Persistence_interface_vineyard.h \ + * gudhi/Multi_persistence/Persistence_interface_vineyard.h + * @ingroup multi_persistence + * + * @brief Interface respecting the @ref PersistenceAlgorithm concept to use @ref Slicer with the homology, vineyard + * and representative cycle algorithms implemented in @ref Gudhi::persistence_matrix::Matrix. + * + * @tparam VineyardOptions Options respecting the @ref VineyardOptions concept. + */ +template +class Persistence_interface_vineyard +{ + public: + using Options = VineyardOptions; + using Vineyard = Gudhi::vineyard::Vineyard_base; /**< Vineyard computation base */ + using Dimension = typename Vineyard::Dimension; /**< Dimension type */ + using Index = typename Vineyard::Index; /**< Index type */ + using Bar = typename Vineyard::Bar; /**< Bar type */ + using Cycle = typename Vineyard::Cycle; /**< Cycle type */ + using Map = typename Vineyard::Permutation; /**< Permutation map for filtration order. */ + template + using As_type = Persistence_interface_vineyard; /**< This type. */ + + static constexpr auto nullDeath = Bar::inf; /**< Death value of the barcode when the bar never died. */ + static constexpr bool is_vine = true; /**< True: update optimization + barcode matching enabled. */ + static constexpr bool has_rep_cycles = true; /**< True: rep. cycle methods enabled. */ + + Persistence_interface_vineyard() = default; + + template + Persistence_interface_vineyard(const Complex& cpx, + const Filtration_range& filtrationValues, + [[maybe_unused]] bool ignoreInf = false) + : vineyard_(cpx.get_boundaries(), cpx.get_dimensions(), filtrationValues) + {} + + Persistence_interface_vineyard(const Persistence_interface_vineyard& other) = default; + + Persistence_interface_vineyard(Persistence_interface_vineyard&& other) noexcept = default; + + ~Persistence_interface_vineyard() = default; + + Persistence_interface_vineyard& operator=(const Persistence_interface_vineyard& other) = default; + + Persistence_interface_vineyard& operator=(Persistence_interface_vineyard&& other) noexcept = default; + + // TODO: swap? + + /** + * @brief The `ignoreInf` argument is ignored for this class, as it disrupts the matching property of the barcodes + * after update. + */ + template + void initialize(const Complex& cpx, const Filtration_range& filtrationValues, [[maybe_unused]] bool ignoreInf = false) + { + vineyard_.initialize(cpx.get_boundaries(), cpx.get_dimensions(), filtrationValues); + } + + /** + * @brief The `ignoreInf` argument is ignored for this class, as it disrupts the matching property of the barcodes + * after update. + */ + template + void update(const Filtration_range& filtrationValues, [[maybe_unused]] bool ignoreInf = false) + { + GUDHI_CHECK(is_initialized(), std::logic_error("Barcode can not be updated uninitialized.")); + + vineyard_.update(filtrationValues); + } + + [[nodiscard]] bool is_initialized() const { return vineyard_.is_initialized(); } + + const Map& get_current_order() const { return vineyard_.get_current_order(); } + + auto get_barcode() + { + GUDHI_CHECK(is_initialized(), std::logic_error("Barcode can not be computed uninitialized.")); + + return vineyard_.get_current_barcode(); + } + + auto get_all_representative_cycles(bool update = true, Dimension dim = Vineyard::nullDimension) + { + GUDHI_CHECK(is_initialized(), std::logic_error("Representative cycles can not be computed uninitialized.")); + + return vineyard_.get_all_current_representative_cycles(update, dim); + } + + auto get_representative_cycle(Index barcodeIndex, bool update = true) + { + GUDHI_CHECK(is_initialized(), std::logic_error("Representative cycles can not be computed uninitialized.")); + + return vineyard_.get_current_representative_cycle(barcodeIndex, update); + } + + friend std::ostream& operator<<(std::ostream& stream, Persistence_interface_vineyard& pers) + { + stream << pers.vineyard_; + + return stream; + } + + private: + Vineyard vineyard_; +}; + +} // namespace multi_persistence +} // namespace Gudhi + +#endif // MP_PERSISTENCE_INTERFACE_VINEYARD_H_INCLUDED diff --git a/multipers/gudhi/gudhi/Persistence_matrix/Base_matrix.h b/multipers/gudhi/gudhi/Persistence_matrix/Base_matrix.h index 5d9c2bc9..83e6a2d9 100644 --- a/multipers/gudhi/gudhi/Persistence_matrix/Base_matrix.h +++ b/multipers/gudhi/gudhi/Persistence_matrix/Base_matrix.h @@ -2,7 +2,7 @@ * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. * Author(s): Hannah Schreiber * - * Copyright (C) 2022-24 Inria + * Copyright (C) 2022 Inria * * Modification(s): * - YYYY/MM Author: Description of the modification diff --git a/multipers/gudhi/gudhi/Persistence_matrix/Base_matrix_with_column_compression.h b/multipers/gudhi/gudhi/Persistence_matrix/Base_matrix_with_column_compression.h index 5246736e..f1440f26 100644 --- a/multipers/gudhi/gudhi/Persistence_matrix/Base_matrix_with_column_compression.h +++ b/multipers/gudhi/gudhi/Persistence_matrix/Base_matrix_with_column_compression.h @@ -2,7 +2,7 @@ * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. * Author(s): Hannah Schreiber * - * Copyright (C) 2022-24 Inria + * Copyright (C) 2022 Inria * * Modification(s): * - YYYY/MM Author: Description of the modification diff --git a/multipers/gudhi/gudhi/Persistence_matrix/Boundary_matrix.h b/multipers/gudhi/gudhi/Persistence_matrix/Boundary_matrix.h index 4d341200..33ea1411 100644 --- a/multipers/gudhi/gudhi/Persistence_matrix/Boundary_matrix.h +++ b/multipers/gudhi/gudhi/Persistence_matrix/Boundary_matrix.h @@ -2,7 +2,7 @@ * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. * Author(s): Hannah Schreiber * - * Copyright (C) 2022-24 Inria + * Copyright (C) 2022 Inria * * Modification(s): * - YYYY/MM Author: Description of the modification diff --git a/multipers/gudhi/gudhi/Persistence_matrix/Chain_matrix.h b/multipers/gudhi/gudhi/Persistence_matrix/Chain_matrix.h index 162bcdb0..5f451046 100644 --- a/multipers/gudhi/gudhi/Persistence_matrix/Chain_matrix.h +++ b/multipers/gudhi/gudhi/Persistence_matrix/Chain_matrix.h @@ -2,7 +2,7 @@ * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. * Author(s): Hannah Schreiber * - * Copyright (C) 2022-24 Inria + * Copyright (C) 2022 Inria * * Modification(s): * - YYYY/MM Author: Description of the modification diff --git a/multipers/gudhi/gudhi/Persistence_matrix/Id_to_index_overlay.h b/multipers/gudhi/gudhi/Persistence_matrix/Id_to_index_overlay.h index ac378193..4f31ba3c 100644 --- a/multipers/gudhi/gudhi/Persistence_matrix/Id_to_index_overlay.h +++ b/multipers/gudhi/gudhi/Persistence_matrix/Id_to_index_overlay.h @@ -2,7 +2,7 @@ * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. * Author(s): Hannah Schreiber * - * Copyright (C) 2022-24 Inria + * Copyright (C) 2022 Inria * * Modification(s): * - YYYY/MM Author: Description of the modification @@ -603,13 +603,13 @@ class Id_to_index_overlay /** * @brief Only available if @ref PersistenceMatrixOptions::can_retrieve_representative_cycles is true. Pre-computes * the representative cycles of the current state of the filtration represented by the matrix. It needs to be called - * before calling @ref get_representative_cycles if the matrix was modified since last call. Otherwise the old cycles - * will be returned. + * before calling @ref get_all_representative_cycles if the matrix was modified since last call. Otherwise the old + * cycles will be returned. * * @param dim If different from default value, only the cycles of the given dimension are updated. * All others are erased. */ - void update_representative_cycles(Dimension dim = Master_matrix::template get_null_value()); + void update_all_representative_cycles(Dimension dim = Master_matrix::template get_null_value()); /** * @brief Only available if @ref PersistenceMatrixOptions::can_retrieve_representative_cycles is true. Pre-computes * the representative cycle in the current matrix state of the given bar. It needs to be called @@ -625,7 +625,7 @@ class Id_to_index_overlay * * @return A const reference to the vector of representative cycles. */ - const std::vector& get_representative_cycles(); + const std::vector& get_all_representative_cycles() const; /** * @brief Only available if @ref PersistenceMatrixOptions::can_retrieve_representative_cycles is true. * Returns the cycle representing the given bar. @@ -633,7 +633,7 @@ class Id_to_index_overlay * @param bar A bar from the current barcode. * @return A const reference to the cycle representing @p bar. */ - const Cycle& get_representative_cycle(const Bar& bar); + const Cycle& get_representative_cycle(const Bar& bar) const; private: using Dictionary = typename Master_matrix::template Dictionary; @@ -1030,9 +1030,9 @@ Id_to_index_overlay::get_current_barcode() } template -inline void Id_to_index_overlay::update_representative_cycles(Dimension dim) +inline void Id_to_index_overlay::update_all_representative_cycles(Dimension dim) { - matrix_.update_representative_cycles(dim); + matrix_.update_all_representative_cycles(dim); } template @@ -1043,14 +1043,14 @@ inline void Id_to_index_overlay::update_repres template inline const std::vector::Cycle>& -Id_to_index_overlay::get_representative_cycles() +Id_to_index_overlay::get_all_representative_cycles() const { - return matrix_.get_representative_cycles(); + return matrix_.get_all_representative_cycles(); } template inline const typename Id_to_index_overlay::Cycle& -Id_to_index_overlay::get_representative_cycle(const Bar& bar) +Id_to_index_overlay::get_representative_cycle(const Bar& bar) const { return matrix_.get_representative_cycle(bar); } diff --git a/multipers/gudhi/gudhi/Persistence_matrix/Position_to_index_overlay.h b/multipers/gudhi/gudhi/Persistence_matrix/Position_to_index_overlay.h index 39373d16..220f7792 100644 --- a/multipers/gudhi/gudhi/Persistence_matrix/Position_to_index_overlay.h +++ b/multipers/gudhi/gudhi/Persistence_matrix/Position_to_index_overlay.h @@ -2,7 +2,7 @@ * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. * Author(s): Hannah Schreiber * - * Copyright (C) 2022-24 Inria + * Copyright (C) 2022 Inria * * Modification(s): * - YYYY/MM Author: Description of the modification @@ -477,13 +477,13 @@ class Position_to_index_overlay /** * @brief Only available if @ref PersistenceMatrixOptions::can_retrieve_representative_cycles is true. Pre-computes * the representative cycles of the current state of the filtration represented by the matrix. It needs to be called - * before calling @ref get_representative_cycles if the matrix was modified since last call. Otherwise the old cycles - * will be returned. + * before calling @ref get_all_representative_cycles if the matrix was modified since last call. Otherwise the old + * cycles will be returned. * * @param dim If different from default value, only the cycles of the given dimension are updated. * All others are erased. */ - void update_representative_cycles(Dimension dim = Master_matrix::template get_null_value()); + void update_all_representative_cycles(Dimension dim = Master_matrix::template get_null_value()); /** * @brief Only available if @ref PersistenceMatrixOptions::can_retrieve_representative_cycles is true. Pre-computes * the representative cycle in the current matrix state of the given bar. It needs to be called @@ -499,7 +499,7 @@ class Position_to_index_overlay * * @return A const reference to the vector of representative cycles. */ - const std::vector& get_representative_cycles(); + const std::vector& get_all_representative_cycles() const; /** * @brief Only available if @ref PersistenceMatrixOptions::can_retrieve_representative_cycles is true. * Returns the cycle representing the given bar. @@ -507,7 +507,7 @@ class Position_to_index_overlay * @param bar A bar from the current barcode. * @return A const reference to the cycle representing @p bar. */ - const Cycle& get_representative_cycle(const Bar& bar); + const Cycle& get_representative_cycle(const Bar& bar) const; /** * @brief Only available if @ref PersistenceMatrixOptions::has_vine_update is true. @@ -814,9 +814,9 @@ Position_to_index_overlay::get_current_barcode } template -inline void Position_to_index_overlay::update_representative_cycles(Dimension dim) +inline void Position_to_index_overlay::update_all_representative_cycles(Dimension dim) { - matrix_.update_representative_cycles(dim); + matrix_.update_all_representative_cycles(dim); } template @@ -827,14 +827,14 @@ inline void Position_to_index_overlay::update_ template inline const std::vector::Cycle>& -Position_to_index_overlay::get_representative_cycles() +Position_to_index_overlay::get_all_representative_cycles() const { - return matrix_.get_representative_cycles(); + return matrix_.get_all_representative_cycles(); } template inline const typename Position_to_index_overlay::Cycle& -Position_to_index_overlay::get_representative_cycle(const Bar& bar) +Position_to_index_overlay::get_representative_cycle(const Bar& bar) const { return matrix_.get_representative_cycle(bar); } diff --git a/multipers/gudhi/gudhi/Persistence_matrix/RU_matrix.h b/multipers/gudhi/gudhi/Persistence_matrix/RU_matrix.h index a7774934..34dd1a60 100644 --- a/multipers/gudhi/gudhi/Persistence_matrix/RU_matrix.h +++ b/multipers/gudhi/gudhi/Persistence_matrix/RU_matrix.h @@ -2,7 +2,7 @@ * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. * Author(s): Hannah Schreiber * - * Copyright (C) 2022-24 Inria + * Copyright (C) 2022 Inria * * Modification(s): * - YYYY/MM Author: Description of the modification diff --git a/multipers/gudhi/gudhi/Persistence_matrix/base_pairing.h b/multipers/gudhi/gudhi/Persistence_matrix/base_pairing.h index e3c5f711..15fabbe5 100644 --- a/multipers/gudhi/gudhi/Persistence_matrix/base_pairing.h +++ b/multipers/gudhi/gudhi/Persistence_matrix/base_pairing.h @@ -2,7 +2,7 @@ * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. * Author(s): Hannah Schreiber * - * Copyright (C) 2022-24 Inria + * Copyright (C) 2022 Inria * * Modification(s): * - YYYY/MM Author: Description of the modification diff --git a/multipers/gudhi/gudhi/Persistence_matrix/base_swap.h b/multipers/gudhi/gudhi/Persistence_matrix/base_swap.h index 66297498..2c8f0c87 100644 --- a/multipers/gudhi/gudhi/Persistence_matrix/base_swap.h +++ b/multipers/gudhi/gudhi/Persistence_matrix/base_swap.h @@ -2,7 +2,7 @@ * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. * Author(s): Hannah Schreiber * - * Copyright (C) 2022-23 Inria + * Copyright (C) 2022 Inria * * Modification(s): * - YYYY/MM Author: Description of the modification diff --git a/multipers/gudhi/gudhi/Persistence_matrix/chain_pairing.h b/multipers/gudhi/gudhi/Persistence_matrix/chain_pairing.h index 470cce2a..e292be89 100644 --- a/multipers/gudhi/gudhi/Persistence_matrix/chain_pairing.h +++ b/multipers/gudhi/gudhi/Persistence_matrix/chain_pairing.h @@ -2,7 +2,7 @@ * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. * Author(s): Hannah Schreiber * - * Copyright (C) 2022-24 Inria + * Copyright (C) 2022 Inria * * Modification(s): * - YYYY/MM Author: Description of the modification diff --git a/multipers/gudhi/gudhi/Persistence_matrix/chain_rep_cycles.h b/multipers/gudhi/gudhi/Persistence_matrix/chain_rep_cycles.h index 7fa1547e..594e1e25 100644 --- a/multipers/gudhi/gudhi/Persistence_matrix/chain_rep_cycles.h +++ b/multipers/gudhi/gudhi/Persistence_matrix/chain_rep_cycles.h @@ -2,7 +2,7 @@ * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. * Author(s): Hannah Schreiber * - * Copyright (C) 2022-24 Inria + * Copyright (C) 2022 Inria * * Modification(s): * - YYYY/MM Author: Description of the modification @@ -24,6 +24,7 @@ #include #endif +#include #include namespace Gudhi { @@ -72,14 +73,14 @@ class Chain_representative_cycles * @param dim If different from default value, only the cycles of the given dimension are updated. * All others are erased. */ - void update_representative_cycles(Dimension dim = Master_matrix::template get_null_value()); + void update_all_representative_cycles(Dimension dim = Master_matrix::template get_null_value()); /** * @brief Computes the current representative cycle of the given bar. All other cycles already computed are left * untouched (and therefore they could be unvalid for the current matrix). * * @note For chain matrices with enabled vine swaps, this method will only be more efficient than - * @ref update_representative_cycles if not called for too many bars. + * @ref update_all_representative_cycles if not called for too many bars. * * @param bar Bar corresponding to the wanted representative cycle. */ @@ -87,20 +88,20 @@ class Chain_representative_cycles /** * @brief Returns the current representative cycles. If the matrix was modified since the last call, - * @ref update_representative_cycles has to be called to update the returned cycles. + * @ref update_all_representative_cycles has to be called to update the returned cycles. * * @return A const reference to a vector of @ref Matrix::Cycle containing all representative cycles. */ - const std::vector& get_representative_cycles(); + const std::vector& get_all_representative_cycles() const; /** * @brief Returns the representative cycle corresponding to the given bar. - * If the matrix was modified since the last call, @ref update_representative_cycles or + * If the matrix was modified since the last call, @ref update_all_representative_cycles or * @ref update_representative_cycle has to be called to update the returned cycle. * * @param bar Bar corresponding to the wanted representative cycle. * @return A const reference to the representative cycle. */ - const Cycle& get_representative_cycle(const Bar& bar); + const Cycle& get_representative_cycle(const Bar& bar) const; /** * @brief Swap operator. @@ -122,12 +123,11 @@ class Chain_representative_cycles // access to inheriting Chain_matrix class constexpr Master_chain_matrix* _matrix() { return static_cast(this); } - constexpr const Master_chain_matrix* _matrix() const { return static_cast(this); } }; template -inline void Chain_representative_cycles::update_representative_cycles(Dimension dim) +inline void Chain_representative_cycles::update_all_representative_cycles(Dimension dim) { Index nberColumns = _matrix()->get_number_of_columns(); auto get_position = [&](ID_index pivot) { @@ -205,31 +205,31 @@ inline void Chain_representative_cycles::update_representative_cy representativeCycles_.resize(representativeCycles_.size() + 1); } + ID_index pivot = bar.birth; + if constexpr (Master_matrix::Option_list::has_vine_update) { - for (ID_index i = 0; i < nberColumns; i++) { - if (get_position(i) == bar.birth) { - auto& col = _matrix()->get_column(_matrix()->get_column_with_pivot(i)); - representativeCycles_[birthToCycle_[bar.birth]] = - Master_matrix::build_cycle_from_range(col.get_non_zero_content_range()); - } - } - } else { - auto& col = _matrix()->get_column(_matrix()->get_column_with_pivot(bar.birth)); - representativeCycles_[birthToCycle_[bar.birth]] = - Master_matrix::build_cycle_from_range(col.get_non_zero_content_range()); + ID_index i = 0; + while (i < nberColumns && get_position(i) != bar.birth) ++i; + GUDHI_CHECK(i < nberColumns, std::invalid_argument("Given birth value is not valid.")); + pivot = i; } + + // if bar.birth does not correspond to a valid birth time, get_column_with_pivot will throw an out of range. + auto& col = _matrix()->get_column(_matrix()->get_column_with_pivot(pivot)); + representativeCycles_[birthToCycle_[bar.birth]] = + Master_matrix::build_cycle_from_range(col.get_non_zero_content_range()); } template inline const std::vector::Cycle>& -Chain_representative_cycles::get_representative_cycles() +Chain_representative_cycles::get_all_representative_cycles() const { return representativeCycles_; } template inline const typename Chain_representative_cycles::Cycle& -Chain_representative_cycles::get_representative_cycle(const Bar& bar) +Chain_representative_cycles::get_representative_cycle(const Bar& bar) const { return representativeCycles_[birthToCycle_[bar.birth]]; } diff --git a/multipers/gudhi/gudhi/Persistence_matrix/chain_vine_swap.h b/multipers/gudhi/gudhi/Persistence_matrix/chain_vine_swap.h index 298216d3..b74dc6b7 100644 --- a/multipers/gudhi/gudhi/Persistence_matrix/chain_vine_swap.h +++ b/multipers/gudhi/gudhi/Persistence_matrix/chain_vine_swap.h @@ -2,7 +2,7 @@ * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. * Author(s): Hannah Schreiber * - * Copyright (C) 2022-24 Inria + * Copyright (C) 2022 Inria * * Modification(s): * - YYYY/MM Author: Description of the modification diff --git a/multipers/gudhi/gudhi/Persistence_matrix/columns/chain_column_extra_properties.h b/multipers/gudhi/gudhi/Persistence_matrix/columns/chain_column_extra_properties.h index 7d18ac73..92477bcf 100644 --- a/multipers/gudhi/gudhi/Persistence_matrix/columns/chain_column_extra_properties.h +++ b/multipers/gudhi/gudhi/Persistence_matrix/columns/chain_column_extra_properties.h @@ -2,7 +2,7 @@ * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. * Author(s): Hannah Schreiber * - * Copyright (C) 2022-24 Inria + * Copyright (C) 2022 Inria * * Modification(s): * - YYYY/MM Author: Description of the modification diff --git a/multipers/gudhi/gudhi/Persistence_matrix/columns/column_dimension_holder.h b/multipers/gudhi/gudhi/Persistence_matrix/columns/column_dimension_holder.h index 72495687..ef67d7c6 100644 --- a/multipers/gudhi/gudhi/Persistence_matrix/columns/column_dimension_holder.h +++ b/multipers/gudhi/gudhi/Persistence_matrix/columns/column_dimension_holder.h @@ -2,7 +2,7 @@ * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. * Author(s): Hannah Schreiber * - * Copyright (C) 2022-24 Inria + * Copyright (C) 2022 Inria * * Modification(s): * - YYYY/MM Author: Description of the modification diff --git a/multipers/gudhi/gudhi/Persistence_matrix/columns/entry_types.h b/multipers/gudhi/gudhi/Persistence_matrix/columns/entry_types.h index 40f9485a..ed46d85c 100644 --- a/multipers/gudhi/gudhi/Persistence_matrix/columns/entry_types.h +++ b/multipers/gudhi/gudhi/Persistence_matrix/columns/entry_types.h @@ -2,7 +2,7 @@ * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. * Author(s): Hannah Schreiber * - * Copyright (C) 2022-24 Inria + * Copyright (C) 2022 Inria * * Modification(s): * - YYYY/MM Author: Description of the modification diff --git a/multipers/gudhi/gudhi/Persistence_matrix/columns/heap_column.h b/multipers/gudhi/gudhi/Persistence_matrix/columns/heap_column.h index 9533d717..86fbbcc3 100644 --- a/multipers/gudhi/gudhi/Persistence_matrix/columns/heap_column.h +++ b/multipers/gudhi/gudhi/Persistence_matrix/columns/heap_column.h @@ -2,7 +2,7 @@ * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. * Author(s): Hannah Schreiber * - * Copyright (C) 2022-24 Inria + * Copyright (C) 2022 Inria * * Modification(s): * - YYYY/MM Author: Description of the modification @@ -772,7 +772,8 @@ inline void Heap_column::push_back(const Entry& entry) { static_assert(Master_matrix::Option_list::is_of_boundary_type, "`push_back` is not available for Chain matrices."); - GUDHI_CHECK(entry.get_row_index() > get_pivot(), "The new row index has to be higher than the current pivot."); + GUDHI_CHECK(entry.get_row_index() > get_pivot(), + std::invalid_argument("The new row index has to be higher than the current pivot.")); Entry* newEntry = entryPool_->construct(entry.get_row_index()); newEntry->set_element(Master_matrix::get_coefficient_value(entry.get_element(), operators_)); diff --git a/multipers/gudhi/gudhi/Persistence_matrix/columns/intrusive_list_column.h b/multipers/gudhi/gudhi/Persistence_matrix/columns/intrusive_list_column.h index 7ec34360..4ccaacd1 100644 --- a/multipers/gudhi/gudhi/Persistence_matrix/columns/intrusive_list_column.h +++ b/multipers/gudhi/gudhi/Persistence_matrix/columns/intrusive_list_column.h @@ -2,7 +2,7 @@ * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. * Author(s): Hannah Schreiber * - * Copyright (C) 2022-24 Inria + * Copyright (C) 2022 Inria * * Modification(s): * - YYYY/MM Author: Description of the modification @@ -807,7 +807,8 @@ inline void Intrusive_list_column::push_back(const Entry& entry) { static_assert(Master_matrix::Option_list::is_of_boundary_type, "`push_back` is not available for Chain matrices."); - GUDHI_CHECK(entry.get_row_index() > get_pivot(), "The new row index has to be higher than the current pivot."); + GUDHI_CHECK(entry.get_row_index() > get_pivot(), + std::invalid_argument("The new row index has to be higher than the current pivot.")); _insert_entry(column_.end(), entry.get_row_index(), entry.get_element()); } diff --git a/multipers/gudhi/gudhi/Persistence_matrix/columns/intrusive_set_column.h b/multipers/gudhi/gudhi/Persistence_matrix/columns/intrusive_set_column.h index 0ee5ded2..2acf8002 100644 --- a/multipers/gudhi/gudhi/Persistence_matrix/columns/intrusive_set_column.h +++ b/multipers/gudhi/gudhi/Persistence_matrix/columns/intrusive_set_column.h @@ -2,7 +2,7 @@ * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. * Author(s): Hannah Schreiber * - * Copyright (C) 2022-24 Inria + * Copyright (C) 2022 Inria * * Modification(s): * - YYYY/MM Author: Description of the modification @@ -595,8 +595,9 @@ Intrusive_set_column::get_pivot_value() const } else { if (Chain_opt::_get_pivot() == Master_matrix::template get_null_value()) return 0; auto it = column_.find(Entry(Chain_opt::_get_pivot())); - GUDHI_CHECK(it != column_.end(), - "Intrusive_set_column::get_pivot_value - Pivot not found only if the column was misused."); + GUDHI_CHECK( + it != column_.end(), + std::logic_error("Intrusive_set_column::get_pivot_value - Pivot not found only if the column was misused.")); return it->get_element(); } } @@ -798,7 +799,8 @@ inline void Intrusive_set_column::push_back(const Entry& entry) { static_assert(Master_matrix::Option_list::is_of_boundary_type, "`push_back` is not available for Chain matrices."); - GUDHI_CHECK(entry.get_row_index() > get_pivot(), "The new row index has to be higher than the current pivot."); + GUDHI_CHECK(entry.get_row_index() > get_pivot(), + std::invalid_argument("The new row index has to be higher than the current pivot.")); _insert_entry(column_.end(), entry.get_row_index(), entry.get_element()); } diff --git a/multipers/gudhi/gudhi/Persistence_matrix/columns/list_column.h b/multipers/gudhi/gudhi/Persistence_matrix/columns/list_column.h index 9699e151..c0c92021 100644 --- a/multipers/gudhi/gudhi/Persistence_matrix/columns/list_column.h +++ b/multipers/gudhi/gudhi/Persistence_matrix/columns/list_column.h @@ -2,7 +2,7 @@ * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. * Author(s): Hannah Schreiber * - * Copyright (C) 2022-24 Inria + * Copyright (C) 2022 Inria * * Modification(s): * - YYYY/MM Author: Description of the modification @@ -775,7 +775,8 @@ inline void List_column::push_back(const Entry& entry) { static_assert(Master_matrix::Option_list::is_of_boundary_type, "`push_back` is not available for Chain matrices."); - GUDHI_CHECK(entry.get_row_index() > get_pivot(), "The new row index has to be higher than the current pivot."); + GUDHI_CHECK(entry.get_row_index() > get_pivot(), + std::invalid_argument("The new row index has to be higher than the current pivot.")); _insert_entry(column_.end(), entry.get_row_index(), entry.get_element()); } diff --git a/multipers/gudhi/gudhi/Persistence_matrix/columns/naive_vector_column.h b/multipers/gudhi/gudhi/Persistence_matrix/columns/naive_vector_column.h index 225d6b22..01942161 100644 --- a/multipers/gudhi/gudhi/Persistence_matrix/columns/naive_vector_column.h +++ b/multipers/gudhi/gudhi/Persistence_matrix/columns/naive_vector_column.h @@ -2,7 +2,7 @@ * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. * Author(s): Hannah Schreiber * - * Copyright (C) 2022-24 Inria + * Copyright (C) 2022 Inria * * Modification(s): * - YYYY/MM Author: Description of the modification @@ -794,7 +794,8 @@ inline void Naive_vector_column::push_back(const Entry& { static_assert(Master_matrix::Option_list::is_of_boundary_type, "`push_back` is not available for Chain matrices."); - GUDHI_CHECK(entry.get_row_index() > get_pivot(), "The new row index has to be higher than the current pivot."); + GUDHI_CHECK(entry.get_row_index() > get_pivot(), + std::invalid_argument("The new row index has to be higher than the current pivot.")); _insert_entry(column_, entry.get_row_index(), entry.get_element()); } diff --git a/multipers/gudhi/gudhi/Persistence_matrix/columns/row_access.h b/multipers/gudhi/gudhi/Persistence_matrix/columns/row_access.h index f3cd461e..7bf32c3c 100644 --- a/multipers/gudhi/gudhi/Persistence_matrix/columns/row_access.h +++ b/multipers/gudhi/gudhi/Persistence_matrix/columns/row_access.h @@ -2,7 +2,7 @@ * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. * Author(s): Hannah Schreiber * - * Copyright (C) 2022-24 Inria + * Copyright (C) 2022 Inria * * Modification(s): * - YYYY/MM Author: Description of the modification diff --git a/multipers/gudhi/gudhi/Persistence_matrix/columns/set_column.h b/multipers/gudhi/gudhi/Persistence_matrix/columns/set_column.h index c1b53523..ca9c7627 100644 --- a/multipers/gudhi/gudhi/Persistence_matrix/columns/set_column.h +++ b/multipers/gudhi/gudhi/Persistence_matrix/columns/set_column.h @@ -2,7 +2,7 @@ * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. * Author(s): Hannah Schreiber * - * Copyright (C) 2022-24 Inria + * Copyright (C) 2022 Inria * * Modification(s): * - YYYY/MM Author: Description of the modification @@ -764,7 +764,8 @@ inline void Set_column::push_back(const Entry& entry) { static_assert(Master_matrix::Option_list::is_of_boundary_type, "`push_back` is not available for Chain matrices."); - GUDHI_CHECK(entry.get_row_index() > get_pivot(), "The new row index has to be higher than the current pivot."); + GUDHI_CHECK(entry.get_row_index() > get_pivot(), + std::invalid_argument("The new row index has to be higher than the current pivot.")); _insert_entry(column_.end(), entry.get_row_index(), entry.get_element()); } diff --git a/multipers/gudhi/gudhi/Persistence_matrix/columns/unordered_set_column.h b/multipers/gudhi/gudhi/Persistence_matrix/columns/unordered_set_column.h index 16997317..be344206 100644 --- a/multipers/gudhi/gudhi/Persistence_matrix/columns/unordered_set_column.h +++ b/multipers/gudhi/gudhi/Persistence_matrix/columns/unordered_set_column.h @@ -2,7 +2,7 @@ * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. * Author(s): Hannah Schreiber * - * Copyright (C) 2022-24 Inria + * Copyright (C) 2022 Inria * * Modification(s): * - YYYY/MM Author: Description of the modification @@ -767,7 +767,8 @@ inline void Unordered_set_column::push_back(const Entry& entry) { static_assert(Master_matrix::Option_list::is_of_boundary_type, "`push_back` is not available for Chain matrices."); - GUDHI_CHECK(entry.get_row_index() > get_pivot(), "The new row index has to be higher than the current pivot."); + GUDHI_CHECK(entry.get_row_index() > get_pivot(), + std::invalid_argument("The new row index has to be higher than the current pivot.")); _insert_entry(entry.get_row_index(), entry.get_element()); } diff --git a/multipers/gudhi/gudhi/Persistence_matrix/columns/vector_column.h b/multipers/gudhi/gudhi/Persistence_matrix/columns/vector_column.h index 8db84f61..888c348d 100644 --- a/multipers/gudhi/gudhi/Persistence_matrix/columns/vector_column.h +++ b/multipers/gudhi/gudhi/Persistence_matrix/columns/vector_column.h @@ -2,7 +2,7 @@ * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. * Author(s): Hannah Schreiber * - * Copyright (C) 2022-24 Inria + * Copyright (C) 2022 Inria * * Modification(s): * - YYYY/MM Author: Description of the modification @@ -865,7 +865,8 @@ inline void Vector_column::push_back(const Entry& entry) { static_assert(Master_matrix::Option_list::is_of_boundary_type, "`push_back` is not available for Chain matrices."); - GUDHI_CHECK(entry.get_row_index() > get_pivot(), "The new row index has to be higher than the current pivot."); + GUDHI_CHECK(entry.get_row_index() > get_pivot(), + std::invalid_argument("The new row index has to be higher than the current pivot.")); _insert_entry(column_, entry.get_row_index(), entry.get_element()); } diff --git a/multipers/gudhi/gudhi/Persistence_matrix/index_mapper.h b/multipers/gudhi/gudhi/Persistence_matrix/index_mapper.h index e3426fea..98128c48 100644 --- a/multipers/gudhi/gudhi/Persistence_matrix/index_mapper.h +++ b/multipers/gudhi/gudhi/Persistence_matrix/index_mapper.h @@ -2,7 +2,7 @@ * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. * Author(s): Hannah Schreiber * - * Copyright (C) 2022-24 Inria + * Copyright (C) 2022 Inria * * Modification(s): * - YYYY/MM Author: Description of the modification diff --git a/multipers/gudhi/gudhi/Persistence_matrix/matrix_dimension_holders.h b/multipers/gudhi/gudhi/Persistence_matrix/matrix_dimension_holders.h index c1720ff2..8210550e 100644 --- a/multipers/gudhi/gudhi/Persistence_matrix/matrix_dimension_holders.h +++ b/multipers/gudhi/gudhi/Persistence_matrix/matrix_dimension_holders.h @@ -2,7 +2,7 @@ * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. * Author(s): Hannah Schreiber * - * Copyright (C) 2022-24 Inria + * Copyright (C) 2022 Inria * * Modification(s): * - YYYY/MM Author: Description of the modification diff --git a/multipers/gudhi/gudhi/Persistence_matrix/matrix_row_access.h b/multipers/gudhi/gudhi/Persistence_matrix/matrix_row_access.h index c1f5be0c..db3c54f2 100644 --- a/multipers/gudhi/gudhi/Persistence_matrix/matrix_row_access.h +++ b/multipers/gudhi/gudhi/Persistence_matrix/matrix_row_access.h @@ -2,7 +2,7 @@ * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. * Author(s): Hannah Schreiber * - * Copyright (C) 2022-24 Inria + * Copyright (C) 2022 Inria * * Modification(s): * - YYYY/MM Author: Description of the modification diff --git a/multipers/gudhi/gudhi/Persistence_matrix/ru_pairing.h b/multipers/gudhi/gudhi/Persistence_matrix/ru_pairing.h index f1f3c1e7..4042f796 100644 --- a/multipers/gudhi/gudhi/Persistence_matrix/ru_pairing.h +++ b/multipers/gudhi/gudhi/Persistence_matrix/ru_pairing.h @@ -2,7 +2,7 @@ * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. * Author(s): Hannah Schreiber * - * Copyright (C) 2022-24 Inria + * Copyright (C) 2022 Inria * * Modification(s): * - YYYY/MM Author: Description of the modification diff --git a/multipers/gudhi/gudhi/Persistence_matrix/ru_rep_cycles.h b/multipers/gudhi/gudhi/Persistence_matrix/ru_rep_cycles.h index b074de3f..e72b3e0f 100644 --- a/multipers/gudhi/gudhi/Persistence_matrix/ru_rep_cycles.h +++ b/multipers/gudhi/gudhi/Persistence_matrix/ru_rep_cycles.h @@ -2,7 +2,7 @@ * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. * Author(s): Hannah Schreiber * - * Copyright (C) 2022-24 Inria + * Copyright (C) 2022 Inria * * Modification(s): * - YYYY/MM Author: Description of the modification @@ -72,7 +72,7 @@ class RU_representative_cycles * @param dim If different from default value, only the cycles of the given dimension are updated. * All others are erased. */ - void update_representative_cycles(Dimension dim = Master_matrix::template get_null_value()); + void update_all_representative_cycles(Dimension dim = Master_matrix::template get_null_value()); /** * @brief Computes the current representative cycle of the given bar. All other cycles already computed are left @@ -84,20 +84,20 @@ class RU_representative_cycles /** * @brief Returns the current representative cycles. If the matrix was modified since the last call, - * @ref update_representative_cycles has to be called to update the returned cycles. + * @ref update_all_representative_cycles has to be called to update the returned cycles. * * @return A const reference to a vector of @ref Matrix::Cycle containing all representative cycles. */ - const std::vector& get_representative_cycles(); + const std::vector& get_all_representative_cycles() const; /** * @brief Returns the representative cycle corresponding to the given bar. - * If the matrix was modified since the last call, @ref update_representative_cycles or + * If the matrix was modified since the last call, @ref update_all_representative_cycles or * @ref update_representative_cycle has to be called to update the returned cycle. * * @param bar Bar corresponding to the wanted representative cycle. * @return A const reference to the representative cycle. */ - const Cycle& get_representative_cycle(const Bar& bar); + const Cycle& get_representative_cycle(const Bar& bar) const; /** * @brief Swap operator. @@ -128,7 +128,7 @@ class RU_representative_cycles }; template -inline void RU_representative_cycles::update_representative_cycles(Dimension dim) +inline void RU_representative_cycles::update_all_representative_cycles(Dimension dim) { Index nberColumns = _matrix()->reducedMatrixR_.get_number_of_columns(); Index nullValue = Master_matrix::template get_null_value(); @@ -195,14 +195,14 @@ inline void RU_representative_cycles::update_representative_cycle template inline const std::vector::Cycle>& -RU_representative_cycles::get_representative_cycles() +RU_representative_cycles::get_all_representative_cycles() const { return representativeCycles_; } template inline const typename RU_representative_cycles::Cycle& -RU_representative_cycles::get_representative_cycle(const Bar& bar) +RU_representative_cycles::get_representative_cycle(const Bar& bar) const { return representativeCycles_[birthToCycle_[bar.birth]]; } diff --git a/multipers/gudhi/gudhi/Persistence_matrix/ru_vine_swap.h b/multipers/gudhi/gudhi/Persistence_matrix/ru_vine_swap.h index 08b2700d..ee532e7f 100644 --- a/multipers/gudhi/gudhi/Persistence_matrix/ru_vine_swap.h +++ b/multipers/gudhi/gudhi/Persistence_matrix/ru_vine_swap.h @@ -2,7 +2,7 @@ * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. * Author(s): Hannah Schreiber * - * Copyright (C) 2022-24 Inria + * Copyright (C) 2022 Inria * * Modification(s): * - YYYY/MM Author: Description of the modification diff --git a/multipers/gudhi/gudhi/Projective_cover_kernel.h b/multipers/gudhi/gudhi/Projective_cover_kernel.h index aff972d7..e1bad3c2 100644 --- a/multipers/gudhi/gudhi/Projective_cover_kernel.h +++ b/multipers/gudhi/gudhi/Projective_cover_kernel.h @@ -36,11 +36,12 @@ namespace multi_persistence { * * @brief TODO (what it is + mention that it only works for 2 parameters) * - * @tparam MultiFiltrationValue Filtration value class respecting the @ref MultiFiltrationValue concept. + * @tparam Complex @ref Multi_parameter_filtered_complex with desired template parameters for the resulting projective + * cover kernel complex. * @tparam columnType Column type to use for the matrix used internally. Default value: * @ref Gudhi::persistence_matrix::Column_types::NAIVE_VECTOR "NAIVE_VECTOR". */ -template class Projective_cover_kernel { @@ -49,18 +50,18 @@ class Projective_cover_kernel * @brief Options for matrix type. */ struct Matrix_options : Gudhi::persistence_matrix::Default_options { - using Index = std::uint32_t; + using Index = typename Complex::Index; + using Dimension = typename Complex::Dimension; }; public: - using Filtration_value = MultiFiltrationValue; /**< Filtration value type. */ - using Complex = Multi_parameter_filtered_complex; /**< Complex data type. */ + using Matrix = Gudhi::persistence_matrix::Matrix; /**< Matrix type. */ + using Filtration_value = typename Complex::Filtration_value; /**< Filtration value type. */ using Index = typename Complex::Index; /**< Index type. */ using Dimension = typename Complex::Dimension; /**< Dimension type. */ using Filtration_value_container = typename Complex::Filtration_value_container; /**< Filt. value container type. */ using Boundary_container = typename Complex::Boundary_container; /**< Boundary container type. */ using Dimension_container = typename Complex::Dimension_container; /**< Dimension container type. */ - using Matrix = Gudhi::persistence_matrix::Matrix; /**< Matrix type. */ // TODO: this only works for 2 parameter modules. Optimize w.r.t. this. /** diff --git a/multipers/gudhi/gudhi/Slicer.h b/multipers/gudhi/gudhi/Slicer.h index 6b16e74b..ddb746b5 100644 --- a/multipers/gudhi/gudhi/Slicer.h +++ b/multipers/gudhi/gudhi/Slicer.h @@ -21,15 +21,19 @@ #include #include #include +#include #include #include -#include //std::iota #include //std::move #include //std::int32_t #include #include // #include //std::stringstream, to remove when to_str gets removed +#ifdef GUDHI_USE_TBB +#include +#endif + #include #include #include @@ -38,7 +42,6 @@ #include #include #include -#include namespace Gudhi { namespace multi_persistence { @@ -63,12 +66,15 @@ template class Slicer { public: - using Persistence = PersistenceAlgorithm; /**< Persistence algorithm type. */ - using Filtration_value = MultiFiltrationValue; /**< Filtration value type. */ - using T = typename Filtration_value::value_type; /**< Numerical filtration value element type. */ - using Complex = Multi_parameter_filtered_complex; /**< Complex type. */ - using Index = typename Complex::Index; /**< Complex index type. */ - using Dimension = typename Complex::Dimension; /**< Dimension type. */ + using Persistence = PersistenceAlgorithm; /**< Persistence algorithm type. */ + using Filtration_value = MultiFiltrationValue; /**< Filtration value type. */ + using T = typename Filtration_value::value_type; /**< Numerical filtration value element type. */ + using Index = typename Persistence::Index; /**< Complex index type. */ + using Dimension = typename Persistence::Dimension; /**< Dimension type. */ + using Complex = Multi_parameter_filtered_complex; /**< Complex type. */ + using Map = typename Persistence::Map; /**< Permutation map for filtration order. */ + using Cycle = std::vector; /**< Cycle type. */ + using Thread_safe = Thread_safe_slicer; /**< Thread safe slicer type. */ template using Bar = Gudhi::persistence_matrix::Persistence_interval; /**< Bar type. */ /** @@ -93,8 +99,6 @@ class Slicer */ template using Multi_dimensional_flat_barcode = std::vector>; - using Cycle = std::vector; /**< Cycle type. */ - using Thread_safe = Thread_safe_slicer; /**< Thread safe slicer type. */ // CONSTRUCTORS @@ -112,7 +116,6 @@ class Slicer Slicer(const Complex& complex) : complex_(complex), slice_(complex.get_number_of_cycle_generators()), - generatorOrder_(complex.get_number_of_cycle_generators()), persistence_() {} @@ -125,7 +128,6 @@ class Slicer Slicer(Complex&& complex) : complex_(std::move(complex)), slice_(complex_.get_number_of_cycle_generators()), - generatorOrder_(complex_.get_number_of_cycle_generators()), persistence_() {} @@ -133,7 +135,7 @@ class Slicer * @brief Copy constructor. Persistence computation initialization is not updated. */ Slicer(const Slicer& other) - : complex_(other.complex_), slice_(other.slice_), generatorOrder_(other.generatorOrder_), persistence_() + : complex_(other.complex_), slice_(other.slice_), persistence_() {} /** @@ -143,7 +145,6 @@ class Slicer Slicer(const Slicer& other) : complex_(other.get_filtered_complex()), slice_(other.get_slice().begin(), other.get_slice().end()), - generatorOrder_(other.get_current_order()), persistence_() {} @@ -153,7 +154,6 @@ class Slicer Slicer(Slicer&& other) noexcept : complex_(std::move(other.complex_)), slice_(std::move(other.slice_)), - generatorOrder_(std::move(other.generatorOrder_)), persistence_() {} @@ -166,8 +166,7 @@ class Slicer { complex_ = other.complex_; slice_ = other.slice_; - generatorOrder_ = other.generatorOrder_; - persistence_.reset(); + persistence_ = Persistence(); return *this; } @@ -180,8 +179,7 @@ class Slicer { complex_ = other.get_filtered_complex(); slice_ = std::vector(other.get_slice().begin(), other.get_slice().end()); - generatorOrder_ = other.get_current_order(); - persistence_.reset(); + persistence_ = Persistence(); return *this; } @@ -193,8 +191,7 @@ class Slicer { complex_ = std::move(other.complex_); slice_ = std::move(other.slice_); - generatorOrder_ = std::move(other.generatorOrder_); - persistence_.reset(); + persistence_ = Persistence(); return *this; } @@ -235,7 +232,7 @@ class Slicer * are not stored in the container. That means that the size can be smaller than what * @ref get_number_of_cycle_generators returns. */ - const std::vector& get_current_order() const { return generatorOrder_; } + const Map& get_current_order() const { return persistence_.get_current_order(); } /** * @brief Returns a const reference to the current slice. It can be initialized or updated with @ref set_slice @@ -393,10 +390,8 @@ class Slicer void prune_above_dimension(int maxDim) { int idx = complex_.prune_above_dimension(maxDim); - generatorOrder_.resize(idx); - generatorOrder_.shrink_to_fit(); slice_.resize(idx); - persistence_.reset(); + persistence_ = Persistence(); } /** @@ -428,49 +423,32 @@ class Slicer * a valid 1-dimensional filtration, the behaviour is undefined. * * @param ignoreInf If true, all cells at infinity filtration values are ignored for the initialization, resulting - * potentially in less storage use and better performance. But note that this can be problematic with the use of - * @ref vineyard_update. Default value: true. + * potentially in less storage use and better performance. But note that this parameter can be potentially ignored + * if the update method of the template parameter @ref PersistenceAlgorithm does not permit this feature. + * Default value: false. */ - void initialize_persistence_computation(const bool ignoreInf = true) + void initialize_persistence_computation(bool ignoreInf = false) { _initialize_persistence_computation(complex_, ignoreInf); } /** - * @brief After the persistence computation was initialized for a slice and the slice changes, this method can - * update everything necessary for the barcode without re-computing everything from scratch (contrary to - * @ref initialize_persistence_computation). Furthermore, it guarantees that the new barcode will "match" the - * precedent one. TODO: explain exactly what it means and how to do the matching. - * The method will have better performance if the complex is ordered by dimension. - * - * Only available if PersistenceAlgorithm::is_vine is true. + * @brief Updates the persistence barcode with the new slice. If @ref PersistenceAlgorithm is + * @ref Persistence_interface_vineyard: will update everything necessary for the barcode without re-computing + * everything from scratch (contrary to @ref initialize_persistence_computation). Furthermore, it guarantees that + * the new barcode will "matches" the precedent one. TODO: explain exactly what it means and how to do the matching. + * For other @ref PersistenceAlgorithm: more or less equivalent to @ref initialize_persistence_computation. * * @pre @ref initialize_persistence_computation has to be called at least once before. * - * @warning If `ignoreInf` was set to true when initializing the persistence computation, any update of the slice has - * to keep at infinity the boundaries which were before, otherwise the behaviour is undefined (it will throw with - * high probability). + * @param ignoreInf If true, all cells at infinity filtration values are ignored in the filtration, resulting + * potentially in less storage use and better performance. But note that this parameter can be potentially ignored + * if the update method of the template parameter `PersistenceAlgorithm` does not permit this feature. + * Default value: false. */ - void vineyard_update() + void update_persistence_computation(bool ignoreInf = false) { - static_assert(Persistence::is_vine, "vineyard_update() not enabled by the chosen PersistenceAlgorithm class."); - - const bool is_ordered_by_dim = complex_.is_ordered_by_dimension(); - // speed up when ordered by dim, to avoid unnecessary swaps - auto dim_condition = [&](int curr) { - if (is_ordered_by_dim) { - return persistence_.get_dimension(curr) == persistence_.get_dimension(curr - 1); - } - return true; - }; - for (Index i = 1; i < generatorOrder_.size(); i++) { - int curr = i; - while (curr > 0 && dim_condition(curr) && slice_[generatorOrder_[curr]] < slice_[generatorOrder_[curr - 1]]) { - persistence_.vine_swap(curr - 1); - std::swap(generatorOrder_[curr - 1], generatorOrder_[curr]); - --curr; - } - } + persistence_.update(slice_, ignoreInf); } /** @@ -544,9 +522,8 @@ class Slicer Index maxIndex = -1; Index maxBirth = std::numeric_limits::max(); T maxLength = 0; - for (Index i = 0; i < barcodeIndices.size(); ++i) { - // barcodeIndices[i] does not work - const auto& bar = barcodeIndices(i); + Index i = 0; + for (const auto& bar : barcodeIndices) { if (bar.dim == dim) { if (bar.death == Persistence::nullDeath) { if (maxBirth > bar.birth) { @@ -562,11 +539,13 @@ class Slicer } } } + ++i; } if (maxIndex == static_cast(-1)) return {}; - return persistence_.get_representative_cycle(maxIndex, update); + auto cycle = persistence_.get_representative_cycle(maxIndex, update); + return {cycle.begin(), cycle.end()}; } // FRIENDS @@ -599,7 +578,7 @@ class Slicer * @brief Builds a new slicer from the given one by projecting its filtration values on a grid. * See @ref coarsen_on_grid with the paramater `coordinate` at true. */ - friend auto build_slicer_coarsen_on_grid(const Slicer& slicer, const std::vector> grid) + friend auto build_slicer_coarsen_on_grid(const Slicer& slicer, const std::vector>& grid) { using return_filtration_value = decltype(std::declval().template as_type()); using return_complex = decltype(build_complex_coarsen_on_grid(slicer.complex_, grid)); @@ -612,7 +591,7 @@ class Slicer */ friend Slicer build_slicer_from_projective_cover_kernel(const Slicer& slicer, Dimension dim) { - Projective_cover_kernel kernel(slicer.complex_, dim); + Projective_cover_kernel kernel(slicer.complex_, dim); return Slicer(kernel.create_complex()); } @@ -652,7 +631,7 @@ class Slicer /** * @brief Outstream operator. */ - friend std::ostream& operator<<(std::ostream& stream, Slicer& slicer) + friend std::ostream& operator<<(std::ostream& stream, const Slicer& slicer) { stream << "-------------------- Slicer \n"; @@ -661,7 +640,7 @@ class Slicer stream << "--- Order \n"; stream << "{"; - for (const auto& idx : slicer.generatorOrder_) stream << idx << ", "; + for (const auto& idx : slicer.get_current_order()) stream << idx << ", "; stream << "}" << '\n'; stream << "--- Current slice filtration\n"; @@ -670,8 +649,8 @@ class Slicer if (!slicer.slice_.empty()) stream << "\b" << "\b"; stream << "}" << '\n'; - stream << "--- PersBackend \n"; - stream << slicer.persistence_; + // stream << "--- PersBackend \n"; + // stream << slicer.persistence_; return stream; } @@ -680,33 +659,34 @@ class Slicer friend Thread_safe; // Thread_safe will use the "_*" methods below instead of "*". // For ThreadSafe version - Slicer(const std::vector& slice, const std::vector& generatorOrder, const Persistence& persistence) - : complex_(), slice_(slice), generatorOrder_(generatorOrder), persistence_(persistence, generatorOrder_) + Slicer(const std::vector& slice, const Persistence& persistence) + : complex_(), slice_(slice), persistence_(persistence) {} - Slicer(std::vector&& slice, std::vector&& generatorOrder, Persistence&& persistence) + Slicer(std::vector&& slice, Persistence&& persistence) : complex_(), slice_(std::move(slice)), - generatorOrder_(std::move(generatorOrder)), - persistence_(std::move(persistence), generatorOrder_) + persistence_(std::move(persistence)) {} template void _push_to(const Complex& complex, const Line& line) { const auto& filtrationValues = complex.get_filtration_values(); - tbb::parallel_for(Index(0), Index(filtrationValues.size()), [&](Index i){ +#ifdef GUDHI_USE_TBB + tbb::parallel_for(Index(0), Index(filtrationValues.size()), [&](Index i) { slice_[i] = line.template compute_forward_intersection(filtrationValues[i]); }); - // for (Index i = 0U; i < filtrationValues.size(); i++) { - // slice_[i] = line.template compute_forward_intersection(filtrationValues[i]); - // } +#else + for (Index i = 0; i < filtrationValues.size(); i++) { + slice_[i] = line.template compute_forward_intersection(filtrationValues[i]); + } +#endif } - void _initialize_persistence_computation(const Complex& complex, const bool ignoreInf = true) + void _initialize_persistence_computation(const Complex& complex, bool ignoreInf = false) { - _initialize_order(complex, ignoreInf); - persistence_.reinitialize(complex, generatorOrder_); + persistence_.initialize(complex, slice_, ignoreInf); } std::vector> _get_representative_cycles(const Complex& complex, bool update = true) @@ -715,7 +695,7 @@ class Slicer "Representative cycles not enabled by the chosen PersistenceAlgorithm class."); const auto& dimensions = complex.get_dimensions(); - auto cycleKeys = persistence_.get_representative_cycles(update); + auto cycleKeys = persistence_.get_all_representative_cycles(update); auto numCycles = cycleKeys.size(); std::vector> out(complex.get_max_dimension() + 1); for (auto& cyclesDim : out) cyclesDim.reserve(numCycles); @@ -730,32 +710,8 @@ class Slicer private: Complex complex_; /**< Complex storing all boundaries, filtration values and dimensions. */ std::vector slice_; /**< Filtration values of the current slice. The indices corresponds to those in complex_. */ - std::vector generatorOrder_; /**< Permutation map from current slice index to complex index. */ Persistence persistence_; /**< Class for persistence computations. */ - void _initialize_order(const Complex& complex, const bool ignoreInf = true) - { - const auto& dimensions = complex.get_dimensions(); - generatorOrder_.resize(complex.get_number_of_cycle_generators()); - std::iota(generatorOrder_.begin(), generatorOrder_.end(), 0); - std::sort(generatorOrder_.begin(), generatorOrder_.end(), [&](Index i, Index j) { - if (ignoreInf) { - if (slice_[i] != Filtration_value::T_inf && slice_[j] == Filtration_value::T_inf) return true; - // all elements at inf are considered equal - if (slice_[i] == Filtration_value::T_inf) return false; - } - if (dimensions[i] > dimensions[j]) return false; - if (dimensions[i] < dimensions[j]) return true; - // if filtration values are equal, we don't care about order, so considered the same object - return slice_[i] < slice_[j]; - }); - if (ignoreInf) { - Index end = generatorOrder_.size(); - while (end > 0 && slice_[generatorOrder_[end - 1]] == Filtration_value::T_inf) --end; - generatorOrder_.resize(end); - } - } - template void _retrieve_interval(const Interval& bar, Dimension& dim, Value& birth, Value& death) { diff --git a/multipers/gudhi/gudhi/Thread_safe_slicer.h b/multipers/gudhi/gudhi/Thread_safe_slicer.h index 115f64e5..d03ad15f 100644 --- a/multipers/gudhi/gudhi/Thread_safe_slicer.h +++ b/multipers/gudhi/gudhi/Thread_safe_slicer.h @@ -78,19 +78,19 @@ class Thread_safe_slicer : private Slicer // CONSTRUCTORS /** - * @brief Constructor. Will store a pointer to the given slicer and copy all persistence related container. + * @brief Constructor. Will store a pointer to the given slicer and copy all persistence related containers. * * @param slicer Original slicer. */ Thread_safe_slicer(const Slicer& slicer) - : Slicer(slicer.get_slice(), slicer.get_current_order(), slicer.get_persistence_algorithm()), slicer_(&slicer) + : Slicer(slicer.get_slice(), slicer.get_persistence_algorithm()), slicer_(&slicer) {} /** * @brief Copy constructor. */ Thread_safe_slicer(const Thread_safe_slicer& slicer) - : Slicer(slicer.get_slice(), slicer.get_current_order(), slicer.get_persistence_algorithm()), + : Slicer(slicer.get_slice(), slicer.get_persistence_algorithm()), slicer_(slicer.slicer_) {} @@ -98,7 +98,7 @@ class Thread_safe_slicer : private Slicer * @brief Move constructor. */ Thread_safe_slicer(Thread_safe_slicer&& slicer) noexcept - : Slicer(std::move(slicer.slice_), std::move(slicer.generatorOrder_), std::move(slicer.persistence_)), + : Slicer(std::move(slicer.slice_), std::move(slicer.persistence_)), slicer_(slicer.slicer_) {} @@ -112,7 +112,6 @@ class Thread_safe_slicer : private Slicer if (this == &other) return *this; Slicer::slice_ = other.slice_; - Slicer::generatorOrder_ = other.generatorOrder_; Slicer::persistence_.reinitialize(slicer_->complex_, Slicer::generatorOrder_); slicer_ = other.slicer_; @@ -270,30 +269,30 @@ class Thread_safe_slicer : private Slicer * a valid 1-dimensional filtration, the behaviour is undefined. * * @param ignoreInf If true, all cells at infinity filtration values are ignored for the initialization, resulting - * potentially in less storage use and better performance. But note that this can be problematic with the use of - * @ref vineyard_update. Default value: true. + * potentially in less storage use and better performance. But note that this parameter can be potentially ignored + * if the update method of the template parameter @ref PersistenceAlgorithm does not permit this feature. + * Default value: false. */ - void initialize_persistence_computation(const bool ignoreInf = true) + void initialize_persistence_computation(bool ignoreInf = false) { Slicer::_initialize_persistence_computation(slicer_->complex_, ignoreInf); } /** - * @brief After the persistence computation was initialized for a slice and the slice changes, this method can - * update everything necessary for the barcode without re-computing everything from scratch (contrary to - * @ref initialize_persistence_computation). Furthermore, it guarantees that the new barcode will "match" the - * precedent one. TODO: explain exactly what it means and how to do the matching. - * The method will have better performance if the complex is ordered by dimension. - * - * Only available if PersistenceAlgorithm::is_vine is true. + * @brief Updates the persistence barcode with the new slice. If @ref PersistenceAlgorithm is + * @ref Persistence_interface_vineyard: will update everything necessary for the barcode without re-computing + * everything from scratch (contrary to @ref initialize_persistence_computation). Furthermore, it guarantees that + * the new barcode will "matches" the precedent one. TODO: explain exactly what it means and how to do the matching. + * For other @ref PersistenceAlgorithm: more or less equivalent to @ref initialize_persistence_computation. * * @pre @ref initialize_persistence_computation has to be called at least once before. * - * @warning If `ignoreInf` was set to true when initializing the persistence computation, any update of the slice has - * to keep at infinity the boundaries which were before, otherwise the behaviour is undefined (it will throw with - * high probability). + * @param ignoreInf If true, all cells at infinity filtration values are ignored in the filtration, resulting + * potentially in less storage use and better performance. But note that this parameter can be potentially ignored + * if the update method of the template parameter `PersistenceAlgorithm` does not permit this feature. + * Default value: false. */ - void vineyard_update() { Slicer::vineyard_update(); } + void update_persistence_computation(bool ignoreInf = false) { Slicer::update_persistence_computation(ignoreInf); } /** * @brief Returns the barcode of the current slice. The barcode format will change depending on the template values. @@ -359,7 +358,7 @@ class Thread_safe_slicer : private Slicer /** * @brief Outstream operator. */ - friend std::ostream& operator<<(std::ostream& stream, Thread_safe_slicer& slicer) + friend std::ostream& operator<<(std::ostream& stream, const Thread_safe_slicer& slicer) { stream << "-------------------- Thread_safe_slicer \n"; @@ -377,8 +376,8 @@ class Thread_safe_slicer : private Slicer stream << "\b" << "\b"; stream << "}" << '\n'; - stream << "--- PersBackend \n"; - stream << slicer.persistence_; + // stream << "--- PersBackend \n"; + // stream << slicer.persistence_; return stream; } diff --git a/multipers/gudhi/gudhi/persistence_matrix_options.h b/multipers/gudhi/gudhi/persistence_matrix_options.h index 3fcadc02..bc1d9f6f 100644 --- a/multipers/gudhi/gudhi/persistence_matrix_options.h +++ b/multipers/gudhi/gudhi/persistence_matrix_options.h @@ -2,7 +2,7 @@ * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. * Author(s): Hannah Schreiber * - * Copyright (C) 2022-24 Inria + * Copyright (C) 2022 Inria * * Modification(s): * - YYYY/MM Author: Description of the modification diff --git a/multipers/gudhi/gudhi/slicer_helpers.h b/multipers/gudhi/gudhi/slicer_helpers.h index 60a4a992..b4379ef1 100644 --- a/multipers/gudhi/gudhi/slicer_helpers.h +++ b/multipers/gudhi/gudhi/slicer_helpers.h @@ -2,7 +2,7 @@ * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. * Author(s): David Loiseaux, Hannah Schreiber * - * Copyright (C) 2023-25 Inria + * Copyright (C) 2023 Inria * * Modification(s): * - YYYY/MM Author: Description of the modification @@ -32,7 +32,6 @@ #include #include #include -#include #ifdef GUDHI_USE_TBB #include @@ -61,6 +60,8 @@ namespace multi_persistence { * * @tparam MultiFiltrationValue Filtration value class respecting the @ref MultiFiltrationValue concept. It will be * used as filtration value type of the new complex. + * @tparam I Index type for the complex. Default value: std::uint32_t. + * @tparam D Dimension type for the complex. Default value: int. * @param inFilePath Path to scc file. * @param isRivetCompatible Set to true if the file is written such that Rivet can read it. See TODO ref. * Default value: false. @@ -71,15 +72,15 @@ namespace multi_persistence { * `shiftDimensions` instead of 0, and if the value is negative, the `abs(shiftDimensions)` smallest dimensions in * the file are ignored and the smallest remaining dimension is interpreted as 0. Default value: 0. */ -template -inline Multi_parameter_filtered_complex build_complex_from_scc_file( +template +inline Multi_parameter_filtered_complex build_complex_from_scc_file( const std::string& inFilePath, bool isRivetCompatible = false, bool isReversed = false, int shiftDimensions = 0) { using Fil = MultiFiltrationValue; - using Complex = Multi_parameter_filtered_complex; + using Complex = Multi_parameter_filtered_complex; using Index = typename Complex::Index; std::string line; @@ -236,6 +237,8 @@ inline Multi_parameter_filtered_complex build_complex_from * boundaries). * * @tparam MultiFiltrationValue Filtration value of the given complex. + * @tparam I Index type of the complex. + * @tparam D Dimension type of the complex. * @param outFilePath Path with file name into which to write. * @param complex Complex to write. Every index appearing in a boundary of the complex has to correspond to an existing * index in the complex @@ -249,9 +252,9 @@ inline Multi_parameter_filtered_complex build_complex_from * @param reverse Set to true if the generators should be written in increasing order of dimension instead of * decreasing. Default value: false. */ -template +template inline void write_complex_to_scc_file(const std::string& outFilePath, - const Multi_parameter_filtered_complex& complex, + const Multi_parameter_filtered_complex& complex, int degree = -1, bool rivetCompatible = false, bool ignoreLastGenerators = false, @@ -298,7 +301,7 @@ inline void write_complex_to_scc_file(const std::string& outFilePath, int minDim = maxDim; const auto& dimensions = complex.get_dimensions(); - std::vector> indicesByDim(maxDim + 1); + std::vector > indicesByDim(maxDim + 1); std::vector shiftedIndices(complex.get_number_of_cycle_generators()); for (std::size_t i = 0; i < complex.get_number_of_cycle_generators(); ++i) { auto dim = dimensions[i]; @@ -397,19 +400,21 @@ inline void write_complex_to_scc_file(const std::string& outFilePath, * * @tparam OneCriticalMultiFiltrationValue Filtration value class respecting the @ref MultiFiltrationValue concept. * It will be used as filtration value type of the new complex. + * @tparam I Index type for the complex. Default value: std::uint32_t. + * @tparam D Dimension type for the complex. Default value: int. * @param vertexValues Bitmap with 1-critical filtration values. Represented as a single vector, the next input * parameter @p shape indicates the shape of the real bitmap. * @param shape Shape of the bitmap. E.g., if @p shape is \f$ {3, 4} \f$, then the bitmap is a \f$ (4 x 3) \f$ grid * with four lines and three columns. The vector @p vertexValues should then contain 12 elements: the three first * elements will be read as the first line, the three next elements as the second line etc. until having 4 lines. */ -template -inline Multi_parameter_filtered_complex build_complex_from_bitmap( +template +inline Multi_parameter_filtered_complex build_complex_from_bitmap( const std::vector& vertexValues, const std::vector& shape) { using Fil = OneCriticalMultiFiltrationValue; - using Complex = Multi_parameter_filtered_complex; + using Complex = Multi_parameter_filtered_complex; using Index = typename Complex::Index; using Bitmap_cubical_complex_base = Gudhi::cubical_complex::Bitmap_cubical_complex_base; using Bitmap_cubical_complex = Gudhi::cubical_complex::Bitmap_cubical_complex; @@ -482,10 +487,12 @@ inline Multi_parameter_filtered_complex build_c * @ref Gudhi::multi_filtration::Multi_parameter_filtration, * @ref Gudhi::multi_filtration::Dynamic_multi_parameter_filtration and * @ref Gudhi::multi_filtration::Degree_rips_bifiltration. + * @tparam I Index type for the complex. Default value: std::uint32_t. + * @tparam D Dimension type for the complex. Default value: int. * @param simplexTree Simplex tree to convert. The key values of the simplex tree will be overwritten. */ -template -inline Multi_parameter_filtered_complex build_complex_from_simplex_tree( +template +inline Multi_parameter_filtered_complex build_complex_from_simplex_tree( Simplex_tree& simplexTree) { // declared here to enable custom `as_type` methods which are not in this namespace. @@ -497,7 +504,7 @@ inline Multi_parameter_filtered_complex build_complex_from static_assert(RangeTraits::is_multi_filtration, "Target filtration value type has to correspond to the MultiFiltrationValue concept."); - using Complex = Multi_parameter_filtered_complex; + using Complex = Multi_parameter_filtered_complex; const unsigned int numberOfSimplices = simplexTree.num_simplices(); @@ -574,7 +581,9 @@ inline Slicer build_slicer_from_scc_file(const std::string& inFilePath, bool isReversed = false, int shiftDimensions = 0) { - auto cpx = build_complex_from_scc_file( + auto cpx = build_complex_from_scc_file( inFilePath, isRivetCompatible, isReversed, shiftDimensions); return Slicer(std::move(cpx)); } @@ -600,7 +609,9 @@ template inline Slicer build_slicer_from_bitmap(const std::vector& vertexValues, const std::vector& shape) { - auto cpx = build_complex_from_bitmap(vertexValues, shape); + auto cpx = + build_complex_from_bitmap( + vertexValues, shape); return Slicer(std::move(cpx)); } @@ -617,7 +628,10 @@ inline Slicer build_slicer_from_bitmap(const std::vector inline Slicer build_slicer_from_simplex_tree(Simplex_tree& simplexTree) { - auto cpx = build_complex_from_simplex_tree(simplexTree); + auto cpx = build_complex_from_simplex_tree(simplexTree); return Slicer(std::move(cpx)); } @@ -625,8 +639,8 @@ inline Slicer build_slicer_from_simplex_tree(Simplex_tree& s * @private */ template -std::vector> -persistence_on_slices_(Slicer& slicer, F&& ini_slicer, unsigned int size, [[maybe_unused]] bool ignoreInf = true) +inline std::vector> +persistence_on_slices_(Slicer& slicer, F&& ini_slicer, unsigned int size, [[maybe_unused]] bool ignoreInf = false) { using Barcode = typename Slicer::template Multi_dimensional_flat_barcode; @@ -640,19 +654,19 @@ persistence_on_slices_(Slicer& slicer, F&& ini_slicer, unsigned int size, [[mayb out[0] = slicer.template get_flat_barcode(); for (auto i = 1U; i < size; ++i) { std::forward(ini_slicer)(slicer, i); - slicer.vineyard_update(); + slicer.update_persistence_computation(); out[i] = slicer.template get_flat_barcode(); } } else { #ifdef GUDHI_USE_TBB - using Index = typename Slicer::Index; tbb::enumerable_thread_specific threadLocals(slicer.weak_copy()); - tbb::parallel_for(static_cast(0), size, [&](const Index& i) { + tbb::parallel_for(static_cast(0), size, [&](const unsigned int& i) { typename Slicer::Thread_safe& s = threadLocals.local(); - // std::forward(ini_slicer)(s, i); - tbb::this_task_arena::isolate([&] { std::forward(ini_slicer)(slicer, i); }); - s.initialize_persistence_computation(ignoreInf); - out[i] = s.template get_flat_barcode(); + tbb::this_task_arena::isolate([&] { + std::forward(ini_slicer)(s, i); // includes another tbb::parallel_for, so needs to be isolated because of s + s.initialize_persistence_computation(ignoreInf); + out[i] = s.template get_flat_barcode(); + }); }); #else for (auto i = 0U; i < size; ++i) { @@ -685,14 +699,14 @@ persistence_on_slices_(Slicer& slicer, F&& ini_slicer, unsigned int size, [[mayb * Can be empty, then the slope is assumed to be 1. * @param ignoreInf If true, all cells at infinity filtration values are ignored when computing, resulting * potentially in less storage use and better performance. But the parameter will be ignored if - * PersistenceAlgorithm::is_vine is true. + * PersistenceAlgorithm::is_vine is true. Default value: false. */ template -std::vector> persistence_on_slices( +inline std::vector> persistence_on_slices( Slicer& slicer, const std::vector>& basePoints, const std::vector>& directions, - bool ignoreInf = true) + bool ignoreInf = false) { GUDHI_CHECK(directions.empty() || directions.size() == basePoints.size(), "There should be as many directions than base points."); @@ -700,13 +714,13 @@ std::vector> persist "There should be as many directions than base points."); std::vector dummy; - auto get_direction = [&](unsigned int i) -> const std::vector& { + auto get_direction = [&](std::size_t i) -> const std::vector& { return directions.empty() ? dummy : directions[i]; }; return persistence_on_slices_( slicer, - [&](auto& s, unsigned int i) { s.push_to(Line(basePoints[i], get_direction(i))); }, + [&](auto& s, std::size_t i) { s.push_to(Line(basePoints[i], get_direction(i))); }, basePoints.size(), ignoreInf); } @@ -725,17 +739,17 @@ std::vector> persist * @param slices Vector of slices. A slice has to has as many elements than cells in the slicer. * @param ignoreInf If true, all cells at infinity filtration values are ignored when computing, resulting * potentially in less storage use and better performance. But the parameter will be ignored if - * PersistenceAlgorithm::is_vine is true. + * PersistenceAlgorithm::is_vine is true. Default value: false. */ template -std::vector> -persistence_on_slices(Slicer& slicer, const std::vector>& slices, bool ignoreInf = true) +inline std::vector> +persistence_on_slices(Slicer& slicer, const std::vector>& slices, bool ignoreInf = false) { GUDHI_CHECK(slices.empty() || slices[0].size() == slicer.get_number_of_cycle_generators(), "There should be as many elements in a slice than cells in the slicer."); return persistence_on_slices_( - slicer, [&](auto& s, unsigned int i) { s.set_slice(slices[i]); }, slices.size(), ignoreInf); + slicer, [&](auto& s, std::size_t i) { s.set_slice(slices[i]); }, slices.size(), ignoreInf); } // Mostly for python @@ -754,18 +768,18 @@ persistence_on_slices(Slicer& slicer, const std::vector>& slices, * @param numberOfSlices Number of slices represented by the pointer. * @param ignoreInf If true, all cells at infinity filtration values are ignored when computing, resulting * potentially in less storage use and better performance. But the parameter will be ignored if - * PersistenceAlgorithm::is_vine is true. + * PersistenceAlgorithm::is_vine is true. Default value: false. */ template >> -std::vector> -persistence_on_slices(Slicer& slicer, T* slices, unsigned int numberOfSlices, bool ignoreInf = true) +inline std::vector> +persistence_on_slices(Slicer& slicer, T* slices, unsigned int numberOfSlices, bool ignoreInf = false) { auto num_gen = slicer.get_number_of_cycle_generators(); auto view = Gudhi::Simple_mdspan(slices, numberOfSlices, num_gen); return persistence_on_slices_( slicer, - [&](auto& s, unsigned int i) { + [&](auto& s, std::size_t i) { T* start = &view(i, 0); auto r = boost::iterator_range(start, start + num_gen); s.set_slice(r); diff --git a/multipers/gudhi/gudhi/vineyard_base.h b/multipers/gudhi/gudhi/vineyard_base.h new file mode 100644 index 00000000..a9b6b9a9 --- /dev/null +++ b/multipers/gudhi/gudhi/vineyard_base.h @@ -0,0 +1,286 @@ +/* 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) 2025 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +/** + * @file vineyard_base.h + * @author Hannah Schreiber + * @brief Contains the implementation of the @ref Gudhi::vineyard::Vineyard_base class. + */ + +#ifndef GUDHI_VINEYARD_BASE_H_ +#define GUDHI_VINEYARD_BASE_H_ + +#include +#include + +#include + +#include +#include + +namespace Gudhi { +namespace vineyard { + +/** + * @ingroup vineyard + * + * @brief Options for the internal matrix of @ref Vineyard_base. + * + * @tparam column_type Column type of the matrix. + */ +template +struct Vineyard_matrix_options : Gudhi::persistence_matrix::Default_options { + using Dimension = D; + using Index = I; + + static const bool is_of_boundary_type = is_RU; + static const Gudhi::persistence_matrix::Column_indexation_types column_indexation_type = + Gudhi::persistence_matrix::Column_indexation_types::POSITION; + static const bool has_column_pairings = true; + static const bool has_vine_update = true; + static const bool can_retrieve_representative_cycles = true; +}; + +/** + * @ingroup vineyard + * + * @brief Default options for @ref Vineyard_base. + */ +struct Default_vineyard_options { + using Index = std::uint32_t; /**< Index type. */ + using Dimension = int; /**< Dimension value type. */ + + static constexpr bool is_RU = false; // TODO: benchmark both + /** + * @brief Column type use by the internal matrix. + */ + static const Gudhi::persistence_matrix::Column_types column_type = + Gudhi::persistence_matrix::Column_types::NAIVE_VECTOR; // TODO: benchmark +}; + +/** + * @class Vineyard_base vineyard_base.h gudhi/vineyard_base.h + * @brief + * @details + * + * @ingroup vineyard + * + * @tparam VineyardOptions Structure following the @ref VineyardOptions concept. Default value: @ref + * Default_vineyard_options. + */ +template +class Vineyard_base +{ + public: + using Index = typename VineyardOptions::Index; /**< Complex index type. */ + using Dimension = typename VineyardOptions::Dimension; /**< Dimension type. */ + using Matrix_options = + Vineyard_matrix_options; + using Matrix = Gudhi::persistence_matrix::Matrix; + using Bar = typename Matrix::Bar; /**< Bar type. */ + using Cycle = typename Matrix::Cycle; /**< Cycle type. */ + using Permutation = std::vector; + + static constexpr Dimension nullDimension = Matrix::template get_null_value(); + + Vineyard_base() = default; + + template + Vineyard_base(const Boundary_range& boundaryMatrix, + const Dimension_range& dimensions, + const Filtration_range& filtrationValues) + : matrix_(boundaryMatrix.size()), order_(boundaryMatrix.size()) + { + // All static_assert in this class are quite useless as Matrix_options is fixed and has those enabled + // I keep them just here for now in case Matrix_options becomes a template argument instead later + static_assert(Matrix_options::has_vine_update, "Underlying matrix has to support vine swaps."); + static_assert(Matrix_options::has_column_pairings, "Underlying matrix has to store barcode."); + static_assert( + Matrix_options::column_indexation_type == Gudhi::persistence_matrix::Column_indexation_types::POSITION || + (Matrix_options::is_of_boundary_type && + Matrix_options::column_indexation_type == Gudhi::persistence_matrix::Column_indexation_types::CONTAINER), + "Matrix has a non supported index scheme."); + + if constexpr (!Matrix_options::is_of_boundary_type) { + idToPos_.emplace(); + idToPos_->resize(order_.size()); + } + + _initialize(boundaryMatrix, dimensions, filtrationValues); + } + + template + void initialize(const Boundary_range& boundaryMatrix, + const Dimension_range& dimensions, + const Filtration_range& filtrationValues) + { + matrix_ = Matrix(boundaryMatrix.size()); + order_.resize(boundaryMatrix.size()); + if constexpr (!Matrix_options::is_of_boundary_type) { + if (!idToPos_) idToPos_.emplace(); + idToPos_->resize(order_.size()); + } + + _initialize(boundaryMatrix, dimensions, filtrationValues); + } + + template + void update(const Filtration_range& filtrationValues) + { + GUDHI_CHECK(filtrationValues.size() == order_.size(), + std::invalid_argument("Filtration value container size is not matching.")); + + for (Index i = 1; i < order_.size(); i++) { + int curr = i; + // speed up when ordered by dim, to avoid unnecessary swaps + while (curr > 0 && matrix_.get_column_dimension(curr) == matrix_.get_column_dimension(curr - 1) && + filtrationValues[order_[curr]] < filtrationValues[order_[curr - 1]]) { + if constexpr (!Matrix_options::is_of_boundary_type) { + auto id1 = matrix_.get_pivot(curr - 1); + auto id2 = matrix_.get_pivot(curr); + std::swap((*idToPos_)[id1], (*idToPos_)[id2]); + } + matrix_.vine_swap(curr - 1); + std::swap(order_[curr - 1], order_[curr]); + --curr; + } + } + } + + [[nodiscard]] bool is_initialized() const { return !order_.empty(); } + + const Permutation& get_current_order() const { return order_; } + + auto get_current_barcode() + { + const auto& barcode = matrix_.get_current_barcode(); + return boost::adaptors::transform(barcode, [&](const Bar& bar) -> Bar { + return Bar(order_[bar.birth], bar.death == Bar::inf ? bar.death : order_[bar.death], bar.dim); + }); + } + + auto get_all_current_representative_cycles(bool update = true, Dimension dim = nullDimension) + { + static_assert(Matrix_options::can_retrieve_representative_cycles, + "Underlying matrix has to support representative cycles."); + + if (update) matrix_.update_all_representative_cycles(dim); + const auto& cycles = matrix_.get_all_representative_cycles(); + return boost::adaptors::transform(cycles, [&](const Cycle& cycle) -> Cycle { + Cycle c(cycle.size()); + for (Index i = 0; i < cycle.size(); ++i) { + if constexpr (Matrix_options::is_of_boundary_type) { + // works for RU because id == pos, but does not work for chain with vine + // we need a id to pos map in that case + c[i] = order_[cycle[i]]; + } else { + c[i] = order_[(*idToPos_)[cycle[i]]]; + } + } + return c; + }); + } + + auto get_current_representative_cycle(Index barcodeIndex, bool update = true) + { + static_assert(Matrix_options::can_retrieve_representative_cycles, + "Underlying matrix has to support representative cycles."); + + const auto& bar = matrix_.get_current_barcode()[barcodeIndex]; + if (update) matrix_.update_representative_cycle(bar); + const auto& cycle = matrix_.get_representative_cycle(bar); + return boost::adaptors::transform(cycle, [&](const Index& i) -> Index { + if constexpr (Matrix_options::is_of_boundary_type) { + // works for RU because id == pos, but does not work for chain with vine + // we need a id to pos map in that case + return order_[i]; + } else { + return order_[(*idToPos_)[i]]; + } + }); + } + + // debug purposes mainly + /** + * @brief Outstream operator. + */ + friend std::ostream& operator<<(std::ostream& stream, Vineyard_base& vyd) + { + stream << "Matrix:\n"; + stream << "[\n"; + for (auto i = 0U; i < vyd.matrix_.get_number_of_columns(); i++) { + stream << "["; + for (const auto& v : vyd.matrix_.get_column(i)) stream << v << ", "; + stream << "]\n"; + } + stream << "]\n"; + stream << "Permutation:\n"; + for (auto v : vyd.order_) { + stream << v << " "; + } + stream << "\n"; + if constexpr (Matrix_options::is_of_boundary_type) { + stream << "ID to position map:\n"; + if (vyd.idToPos_) { + for (auto v : *vyd.idToPos_) { + stream << v << " "; + } + } + stream << "\n"; + } + return stream; + } + + private: + Matrix matrix_; + Permutation order_; + std::optional idToPos_; // TODO: remove if chain does not improve run times in benchmark + + template + void _initialize(const Boundary_range& boundaryMatrix, + const Dimension_range& dimensions, + const Filtration_range& filtrationValues) + { + GUDHI_CHECK(boundaryMatrix.size() == dimensions.size(), + std::invalid_argument("Boundary and dimension range sizes are not matching.")); + GUDHI_CHECK(boundaryMatrix.size() == filtrationValues.size(), + std::invalid_argument("Boundary and filtration value range sizes are not matching.")); + + std::iota(order_.begin(), order_.end(), 0); + std::sort(order_.begin(), order_.end(), [&](Index i, Index j) { + if (dimensions[i] < dimensions[j]) return true; + if (dimensions[i] > dimensions[j]) return false; + return filtrationValues[i] < filtrationValues[j]; + }); + + // simplex IDs need to be increasing in order, so the original ones cannot be used + Permutation orderInv(boundaryMatrix.size()); + Permutation translatedBoundary; + Index id = 0; + for (auto i : order_) { + orderInv[i] = id; // order is assumed to be a valid filtration + translatedBoundary.resize(boundaryMatrix[i].size()); + for (Index j = 0; j < boundaryMatrix[i].size(); ++j) { + translatedBoundary[j] = orderInv[boundaryMatrix[i][j]]; + } + std::sort(translatedBoundary.begin(), translatedBoundary.end()); + matrix_.insert_boundary(id, translatedBoundary, dimensions[i]); + if constexpr (!Matrix_options::is_of_boundary_type) { + (*idToPos_)[id] = id; // before any vine swaps, id == pos + } + ++id; // IDs corresponds to the indices in order_ + } + } +}; + +} // namespace vineyard +} // namespace Gudhi + +#endif // GUDHI_VINEYARD_BASE_H_ diff --git a/multipers/io.pyx b/multipers/io.pyx index e14afa60..26766644 100644 --- a/multipers/io.pyx +++ b/multipers/io.pyx @@ -327,37 +327,37 @@ def _multi_critical_from_slicer( **slicer_kwargs ) with tempfile.TemporaryDirectory(prefix="multipers", delete=clear) as tmpdir: - input_path = os.path.join(tmpdir, "multipers_input.scc") - output_path = os.path.join(tmpdir, "multipers_output.scc") - slicer.to_scc(input_path, degree=0, strip_comments=True) - reduce_arg = "" - if reduce: - if swedish: - reduce_arg += r" --swedish" - if degree is None: - need_split = True - reduce_arg += r" --minpres-all" - else: - reduce_arg += fr" --minpres {slicer.dimension - degree+2}" - verbose_arg = "> /dev/null 2>&1" if not verbose else "" - - command = f"{pathes['multi_critical']} --{algo} {reduce_arg} {input_path} {output_path} {verbose_arg}" - if verbose: - print(command) - os.system(command) - if need_split: - os.system(f"awk \'/scc2020/ {{n++}} {{print > (\"{tmpdir}/multipers_block_\" n \".scc\")}}\' {output_path}") - from glob import glob - import re - files = glob(tmpdir + "/multipers_block_*.scc") - files.sort(key=lambda f: int(re.search(r'\d+', f).group())) - num_degrees=len(files) - ss = tuple(newSlicer()._build_from_scc_file(files[i], shift_dimension=i-1).minpres(i) for i in range(num_degrees)) - return ss - out = newSlicer()._build_from_scc_file(str(output_path), shift_dimension=degree-1 if reduce else -2) - if reduce: - out = out.minpres(degree) - return out + input_path = os.path.join(tmpdir, "multipers_input.scc") + output_path = os.path.join(tmpdir, "multipers_output.scc") + slicer.to_scc(input_path, degree=0, strip_comments=True) + reduce_arg = "" + if reduce: + if swedish: + reduce_arg += r" --swedish" + if degree is None: + need_split = True + reduce_arg += r" --minpres-all" + else: + reduce_arg += fr" --minpres {slicer.dimension - degree+2}" + verbose_arg = "> /dev/null 2>&1" if not verbose else "" + + command = f"{pathes['multi_critical']} --{algo} {reduce_arg} {input_path} {output_path} {verbose_arg}" + if verbose: + print(command) + os.system(command) + if need_split: + os.system(f"awk \'/scc2020/ {{n++}} {{print > (\"{tmpdir}/multipers_block_\" n \".scc\")}}\' {output_path}") + from glob import glob + import re + files = glob(tmpdir + "/multipers_block_*.scc") + files.sort(key=lambda f: int(re.search(r'\d+', f).group())) + num_degrees=len(files) + ss = tuple(newSlicer()._build_from_scc_file(files[i], shift_dimension=i-1).minpres(i) for i in range(num_degrees)) + return ss + out = newSlicer()._build_from_scc_file(str(output_path), shift_dimension=degree-1 if reduce else -2) + if reduce: + out = out.minpres(degree) + return out diff --git a/multipers/multi_parameter_rank_invariant/hilbert_function.h b/multipers/multi_parameter_rank_invariant/hilbert_function.h index b2a632f4..a9df39ca 100644 --- a/multipers/multi_parameter_rank_invariant/hilbert_function.h +++ b/multipers/multi_parameter_rank_invariant/hilbert_function.h @@ -435,14 +435,10 @@ inline void compute_2d_hilbert_surface( using bc_type = typename Gudhi::multi_persistence::Slicer::template Multi_dimensional_flat_barcode<>; - if constexpr (PersBackend::is_vine) { - if (!slicer.persistence_computation_is_initialized()) [[unlikely]] { - slicer.initialize_persistence_computation(); - } else { - slicer.vineyard_update(); - } - } else { + if (!slicer.persistence_computation_is_initialized()) [[unlikely]] { slicer.initialize_persistence_computation(ignore_inf); + } else { + slicer.update_persistence_computation(ignore_inf); } bc_type barcodes = slicer.template get_flat_barcode(); diff --git a/multipers/multi_parameter_rank_invariant/rank_invariant.h b/multipers/multi_parameter_rank_invariant/rank_invariant.h index 7d2874f4..ce05ccf9 100644 --- a/multipers/multi_parameter_rank_invariant/rank_invariant.h +++ b/multipers/multi_parameter_rank_invariant/rank_invariant.h @@ -212,7 +212,7 @@ inline void compute_2d_rank_invariant_of_elbow( // slicer.compute_persistence(degrees_index); slicer.initialize_persistence_computation(); } else { - slicer.vineyard_update(); + slicer.update_persistence_computation(); } } else { slicer.initialize_persistence_computation(ignore_inf); diff --git a/multipers/multiparameter_module_approximation/approximation.h b/multipers/multiparameter_module_approximation/approximation.h index 43642845..249b6b03 100644 --- a/multipers/multiparameter_module_approximation/approximation.h +++ b/multipers/multiparameter_module_approximation/approximation.h @@ -556,7 +556,7 @@ inline void __add_vineyard_trajectory_to_module(Module(), threshold); }; @@ -586,7 +586,7 @@ void _rec_mma(Module &module, _rec_mma(module, basepoint_copy, grid_size, dim_to_iterate - 1, pers_copy, precision, threshold); basepoint[dim_to_iterate] += precision; // current_persistence.push_to(Line(basepoint)); - // current_persistence.vineyard_update(); + // current_persistence.update_persistence_computation(); } } @@ -645,7 +645,7 @@ void _rec_mma2(Module &module, threshold); basepoint[dim_to_iterate] += signs[dim_to_iterate] ? precision : -precision; // current_persistence.push_to(Line(basepoint)); - // current_persistence.vineyard_update(); + // current_persistence.update_persistence_computation(); } } diff --git a/multipers/slicer.pxd.tp b/multipers/slicer.pxd.tp index 0035a9b8..63d95758 100644 --- a/multipers/slicer.pxd.tp +++ b/multipers/slicer.pxd.tp @@ -76,7 +76,6 @@ cdef extern from "Persistence_slices_interface.h" namespace "multipers::tmp_inte uint32_t num_parameters "get_number_of_parameters"() nogil {{if D['IS_VINE']}} vector[uint32_t] get_current_order() nogil - void vineyard_update() nogil vector[vector[vector[uint32_t]]] get_representative_cycles(bool) nogil vector[uint32_t] get_most_persistent_cycle(int, bool) nogil {{endif}} @@ -92,7 +91,9 @@ cdef extern from "Persistence_slices_interface.h" namespace "multipers::tmp_inte int prune_above_dimension(int) except + nogil void coarsen_on_grid_inplace "coarsen_on_grid"(vector[vector[{{D['C_VALUE_TYPE']}}]], bool) nogil void initialize_persistence_computation(bool) except+ nogil - void initialize_persistence_computation() except+ nogil # ignore_inf = true + void initialize_persistence_computation() except+ nogil # ignore_inf = false + void update_persistence_computation(bool) except+ nogil + void update_persistence_computation() except+ nogil # ignore_inf = false Dim_barcode[{{D['C_TEMPLATE_TYPE']}}, {{D['C_VALUE_TYPE']}}] get_barcode "get_flat_barcode"() nogil Dim_barcode[{{D['C_TEMPLATE_TYPE']}}, int] get_barcode_idx "get_flat_barcode"() nogil string slicer_to_str({{D['C_TEMPLATE_TYPE']}}&) nogil #to_str diff --git a/multipers/slicer.pyx.tp b/multipers/slicer.pyx.tp index 8ffb3654..3363f303 100644 --- a/multipers/slicer.pyx.tp +++ b/multipers/slicer.pyx.tp @@ -792,7 +792,7 @@ cdef class {{D['PYTHON_TYPE']}}: Updates the barcode, on a line, using the vineyard algorithm. """ self.push_to_line(basepoint,direction) - self.truc.vineyard_update() + self.truc.update_persistence_computation() return self def get_representative_cycles(self, bool update=True): """