diff --git a/multipers/gudhi/gudhi/Matrix.h b/multipers/gudhi/gudhi/Matrix.h index d917bdd..82905b3 100644 --- a/multipers/gudhi/gudhi/Matrix.h +++ b/multipers/gudhi/gudhi/Matrix.h @@ -1220,7 +1220,7 @@ class Matrix { std::swap(matrix1.colSettings_, matrix2.colSettings_); } - void print(); // for debug + void print(Index startCol = 0, Index endCol = -1, Index startRow = 0, Index endRow = -1); // for debug //TODO: change the behaviour for boundary matrices. /** @@ -1344,21 +1344,23 @@ class Matrix { * returned. */ void update_representative_cycles(); - /** - * @brief Only available if @ref PersistenceMatrixOptions::can_retrieve_representative_cycles is true. - * Returns all representative cycles of the current filtration. - * - * @return A const reference to the vector of representative cycles. - */ - const std::vector& get_representative_cycles(); - /** - * @brief Only available if @ref PersistenceMatrixOptions::can_retrieve_representative_cycles is true. - * Returns the cycle representing the given bar. - * - * @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); + // /** + // * @brief Only available if @ref PersistenceMatrixOptions::can_retrieve_representative_cycles is true. + // * Returns all representative cycles of the current filtration. + // * + // * @return A const reference to the vector of representative cycles. + // */ + // const std::vector& get_representative_cycles(); + // /** + // * @brief Only available if @ref PersistenceMatrixOptions::can_retrieve_representative_cycles is true. + // * Returns the cycle representing the given bar. + // * + // * @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); + std::vector::ID_index> > > + get_representative_cycles_as_borders(bool detailed = false); private: using Underlying_matrix = @@ -1931,9 +1933,9 @@ inline Matrix& Matrix::opera } template -inline void Matrix::print() +inline void Matrix::print(Index startCol, Index endCol, Index startRow, Index endRow) { - return matrix_.print(); + return matrix_.print(startCol, endCol, startRow, endRow); } template @@ -2035,20 +2037,28 @@ inline void Matrix::update_representative_cycles() matrix_.update_representative_cycles(); } -template -inline const std::vector::Cycle>& -Matrix::get_representative_cycles() -{ - static_assert(PersistenceMatrixOptions::can_retrieve_representative_cycles, "This method was not enabled."); - return matrix_.get_representative_cycles(); -} +// template +// inline const std::vector::Cycle>& +// Matrix::get_representative_cycles() +// { +// static_assert(PersistenceMatrixOptions::can_retrieve_representative_cycles, "This method was not enabled."); +// return matrix_.get_representative_cycles(); +// } + +// template +// inline const typename Matrix::Cycle& +// Matrix::get_representative_cycle(const Bar& bar) +// { +// static_assert(PersistenceMatrixOptions::can_retrieve_representative_cycles, "This method was not enabled."); +// return matrix_.get_representative_cycle(bar); +// } template -inline const typename Matrix::Cycle& -Matrix::get_representative_cycle(const Bar& bar) -{ - static_assert(PersistenceMatrixOptions::can_retrieve_representative_cycles, "This method was not enabled."); - return matrix_.get_representative_cycle(bar); +inline std::vector::ID_index> > > +Matrix::get_representative_cycles_as_borders(bool detailed) { + static_assert(PersistenceMatrixOptions::can_retrieve_representative_cycles && PersistenceMatrixOptions::is_z2, + "This method was not enabled."); + return matrix_.get_representative_cycles_as_borders(detailed); } template diff --git a/multipers/gudhi/gudhi/Persistence_matrix/Base_matrix.h b/multipers/gudhi/gudhi/Persistence_matrix/Base_matrix.h index 0630aab..ba8f181 100644 --- a/multipers/gudhi/gudhi/Persistence_matrix/Base_matrix.h +++ b/multipers/gudhi/gudhi/Persistence_matrix/Base_matrix.h @@ -316,7 +316,7 @@ class Base_matrix : public Master_matrix::template Base_swap_option >; @@ -597,27 +597,32 @@ inline Base_matrix& Base_matrix::operator=(const B } template -inline void Base_matrix::print() +inline void Base_matrix::print(Index startCol, Index endCol, Index startRow, Index endRow) { _orderRowsIfNecessary(); + if (endCol == static_cast(-1)) endCol = nextInsertIndex_; + if (endRow == static_cast(-1)) endRow = nextInsertIndex_; std::cout << "Base_matrix:\n"; - for (Index i = 0; i < nextInsertIndex_; ++i) { + for (Index i = startCol; i < endCol && i < nextInsertIndex_; ++i) { const Column& col = matrix_[i]; - for (const auto& e : col.get_content(nextInsertIndex_)) { - if (e == 0u) + auto cont = col.get_content(endRow); + for (Index j = startRow; j < endRow; ++j) { + if (cont[j] == 0u) std::cout << "- "; else - std::cout << e << " "; + std::cout << cont[j] << " "; } std::cout << "\n"; } std::cout << "\n"; if constexpr (Master_matrix::Option_list::has_row_access) { std::cout << "Row Matrix:\n"; - for (Index i = 0; i < nextInsertIndex_; ++i) { + for (Index i = startRow; i < endRow && i < nextInsertIndex_; ++i) { const auto& row = (*RA_opt::rows_)[i]; - for (const auto& entry : row) { - std::cout << entry.get_column_index() << " "; + auto it = row.begin(); + std::advance(it, startCol); + for (; it != row.end() && startCol < endCol; ++it, ++startCol) { + std::cout << it->get_column_index() << " "; } std::cout << "(" << i << ")\n"; } diff --git a/multipers/gudhi/gudhi/Persistence_matrix/Boundary_matrix.h b/multipers/gudhi/gudhi/Persistence_matrix/Boundary_matrix.h index 4ffb3fb..3330d02 100644 --- a/multipers/gudhi/gudhi/Persistence_matrix/Boundary_matrix.h +++ b/multipers/gudhi/gudhi/Persistence_matrix/Boundary_matrix.h @@ -19,6 +19,7 @@ #include #include //print() only +#include #include #include //std::swap, std::move & std::exchange @@ -349,7 +350,7 @@ class Boundary_matrix : public Master_matrix::Matrix_dimension_option, } } - void print(); // for debug + void print(Index startCol = 0, Index endCol = -1, Index startRow = 0, Index endRow = -1); // for debug private: using Dim_opt = typename Master_matrix::Matrix_dimension_option; @@ -709,29 +710,34 @@ inline Boundary_matrix& Boundary_matrix::operator= } template -inline void Boundary_matrix::print() +inline void Boundary_matrix::print(Index startCol, Index endCol, Index startRow, Index endRow) { + if (endCol == static_cast(-1)) endCol = nextInsertIndex_; + if (endRow == static_cast(-1)) endRow = nextInsertIndex_; if constexpr (activeSwapOption) { if (Swap_opt::rowSwapped_) Swap_opt::_orderRows(); } std::cout << "Boundary_matrix:\n"; - for (Index i = 0; i < nextInsertIndex_; ++i) { + for (Index i = startCol; i < endCol && i < nextInsertIndex_; ++i) { Column& col = matrix_[i]; - for (auto e : col.get_content(nextInsertIndex_)) { - if (e == 0u) + auto cont = col.get_content(endRow); + for (Index j = startRow; j < endRow; ++j) { + if (cont[j] == 0u) std::cout << "- "; else - std::cout << e << " "; + std::cout << cont[j] << " "; } std::cout << "\n"; } std::cout << "\n"; if constexpr (Master_matrix::Option_list::has_row_access) { std::cout << "Row Matrix:\n"; - for (ID_index i = 0; i < nextInsertIndex_; ++i) { + for (ID_index i = startRow; i < endRow && i < nextInsertIndex_; ++i) { const auto& row = (*RA_opt::rows_)[i]; - for (const typename Column::Entry& entry : row) { - std::cout << entry.get_column_index() << " "; + auto it = row.begin(); + std::advance(it, startCol); + for (; it != row.end() && startCol < endCol; ++it, ++startCol) { + std::cout << it->get_column_index() << " "; } std::cout << "(" << i << ")\n"; } diff --git a/multipers/gudhi/gudhi/Persistence_matrix/RU_matrix.h b/multipers/gudhi/gudhi/Persistence_matrix/RU_matrix.h index 38ceae6..2b2a336 100644 --- a/multipers/gudhi/gudhi/Persistence_matrix/RU_matrix.h +++ b/multipers/gudhi/gudhi/Persistence_matrix/RU_matrix.h @@ -380,7 +380,7 @@ class RU_matrix : public Master_matrix::RU_pairing_option, std::swap(matrix1.operators_, matrix2.operators_); } - void print(); // for debug + void print(Index startCol = 0, Index endCol = -1, Index startRow = 0, Index endRow = -1); // for debug private: using Swap_opt = typename Master_matrix::RU_vine_swap_option; @@ -538,7 +538,7 @@ inline void RU_matrix::insert_boundary(ID_index cellIndex, if constexpr (Master_matrix::Option_list::has_vine_update) { if (cellIndex != nextEventIndex_) { Swap_opt::_positionToRowIdx().emplace(nextEventIndex_, cellIndex); - if (Master_matrix::Option_list::has_column_pairings) { + if constexpr (Master_matrix::Option_list::has_column_pairings) { Swap_opt::template RU_pairing::idToPosition_.emplace(cellIndex, nextEventIndex_); } } @@ -733,12 +733,12 @@ inline RU_matrix& RU_matrix::operator=(const RU_ma } template -inline void RU_matrix::print() +inline void RU_matrix::print(Index startCol, Index endCol, Index startRow, Index endRow) { std::cout << "R_matrix:\n"; - reducedMatrixR_.print(); + reducedMatrixR_.print(startCol, endCol, startRow, endRow); std::cout << "U_matrix:\n"; - mirrorMatrixU_.print(); + mirrorMatrixU_.print(startCol, endCol, startCol, endCol); } template diff --git a/multipers/gudhi/gudhi/Persistence_matrix/allocators/cell_constructors.h b/multipers/gudhi/gudhi/Persistence_matrix/allocators/cell_constructors.h deleted file mode 100644 index f13d922..0000000 --- a/multipers/gudhi/gudhi/Persistence_matrix/allocators/cell_constructors.h +++ /dev/null @@ -1,139 +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) 2024 Inria - * - * Modification(s): - * - YYYY/MM Author: Description of the modification - */ - -/** - * @file cell_constructors.h - * @author Hannah Schreiber - * @brief Contains different versions of @ref Gudhi::persistence_matrix::Cell factories. - */ - -#ifndef PM_COLUMN_CELL_CONSTRUCTORS_H -#define PM_COLUMN_CELL_CONSTRUCTORS_H - -#include //std::swap - -#include - -namespace Gudhi { -namespace persistence_matrix { - -/** - * @private - * @ingroup persistence_matrix - * - * @brief @ref Cell factory. Constructs and destroys cell pointers with new and delete. - * - * @tparam Cell @ref Cell with the right templates. - */ -template -struct New_cell_constructor -{ - /** - * @brief Default constructor. - */ - New_cell_constructor() {} - - /** - * @brief Constructs a cell with the given cell arguments. - * - * @param u Arguments forwarded to the @ref Cell constructor. - * @return @ref Cell pointer. - */ - template - Cell* construct(U&&... u) const { - return new Cell(std::forward(u)...); - } - - /** - * @brief Destroys the given cell. - * - * @param cell @ref Cell pointer. - */ - void destroy(Cell* cell) const { delete cell; } - - /** - * @brief Swap operator. - */ - friend void swap(New_cell_constructor& col1, New_cell_constructor& col2) {} -}; - -/** - * @private - * @ingroup persistence_matrix - * - * @brief @ref Cell factory. Uses @ref Gudhi::Simple_object_pool, which is based on boost::object_pool, - * to construct and destroy cell pointer. - * - * @tparam Cell @ref Cell with the right templates. - */ -template -struct Pool_cell_constructor -{ - public: - /** - * @brief Default constructor. - * - */ - Pool_cell_constructor() : cellPool_() {} - //TODO: what does happen when the pool is copied? - /** - * @brief Copy constructor. - * - * @param col Factory to copy. - */ - Pool_cell_constructor(const Pool_cell_constructor& col) : cellPool_(col.cellPool_) {} - /** - * @brief Move constructor. - * - * @param col Factory to move. - */ - Pool_cell_constructor(Pool_cell_constructor&& col) : cellPool_(std::move(col.cellPool_)) {} - - /** - * @brief Constructs a cell with the given cell arguments. - * - * @param u Arguments forwarded to the @ref Cell constructor. - * @return @ref Cell pointer. - */ - template - Cell* construct(U&&... u) { - return cellPool_.construct(std::forward(u)...); - } - - /** - * @brief Destroys the given cell. - * - * @param cell @ref Cell pointer. - */ - void destroy(Cell* cell) { cellPool_.destroy(cell); } - - //TODO: Again, what does it mean to copy the pool? - /** - * @brief Assign operator. - */ - Pool_cell_constructor& operator=(const Pool_cell_constructor& other) { - cellPool_ = other.cellPool_; - return *this; - } - /** - * @brief Swap operator. - */ - friend void swap(Pool_cell_constructor& col1, Pool_cell_constructor& col2) { - std::swap(col1.cellPool_, col2.cellPool_); - } - - private: - Simple_object_pool cellPool_; /**< Cell pool. */ -}; - -} // namespace persistence_matrix -} // namespace Gudhi - -#endif // PM_COLUMN_CELL_CONSTRUCTORS_H diff --git a/multipers/gudhi/gudhi/Persistence_matrix/columns/cell_types.h b/multipers/gudhi/gudhi/Persistence_matrix/columns/cell_types.h deleted file mode 100644 index 0bf59de..0000000 --- a/multipers/gudhi/gudhi/Persistence_matrix/columns/cell_types.h +++ /dev/null @@ -1,327 +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) 2022-24 Inria - * - * Modification(s): - * - YYYY/MM Author: Description of the modification - */ - -/** - * @file cell_types.h - * @author Hannah Schreiber - * @brief Contains the @ref Gudhi::persistence_matrix::Cell, @ref Gudhi::persistence_matrix::Cell_column_index and - * @ref Gudhi::persistence_matrix::Cell_field_element classes, as well as the - * @ref Gudhi::persistence_matrix::Dummy_cell_column_index_mixin and - * @ref Gudhi::persistence_matrix::Dummy_cell_field_element_mixin structures. - * Also defines the std::hash method for @ref Gudhi::persistence_matrix::Cell. - */ - -#ifndef PM_MATRIX_CELL_H -#define PM_MATRIX_CELL_H - -#include //std::swap, std::exchange & std::move -#include //std::hash - -namespace Gudhi { -namespace persistence_matrix { - -/** - * @ingroup persistence_matrix - * - * @brief Empty structure. - * Inherited instead of @ref Cell_column_index, when the row access is disabled. - */ -struct Dummy_cell_column_index_mixin -{ - Dummy_cell_column_index_mixin() {} - template - Dummy_cell_column_index_mixin([[maybe_unused]] Index columnIndex) {} -}; - -/** - * @ingroup persistence_matrix - * - * @brief Empty structure. - * Inherited instead of @ref Cell_field_element, when @ref PersistenceMatrixOptions::is_z2 is true. - */ -struct Dummy_cell_field_element_mixin -{ - Dummy_cell_field_element_mixin() {} - template - Dummy_cell_field_element_mixin([[maybe_unused]] Field_element t) {} -}; - -/** - * @ingroup persistence_matrix - * - * @brief Class managing the column index access of a cell. - * - * @tparam Index @ref MatIdx index type. - */ -template -class Cell_column_index -{ - public: - /** - * @brief Default constructor. Sets to the column index to -1. - */ - Cell_column_index() : columnIndex_(-1){}; - /** - * @brief Stores the given column index. - * - * @param columnIndex Column index of the cell. - */ - Cell_column_index(Index columnIndex) : columnIndex_(columnIndex){}; - /** - * @brief Copy constructor. - * - * @param cell Cell to copy. - */ - Cell_column_index(const Cell_column_index& cell) : columnIndex_(cell.columnIndex_){}; - /** - * @brief Move constructor. - * - * @param cell Cell to move. - */ - Cell_column_index(Cell_column_index&& cell) noexcept : columnIndex_(std::exchange(cell.columnIndex_, 0)){}; - - /** - * @brief Returns the @ref MatIdx column index stored in the cell. - * - * @return Column index of the cell. - */ - Index get_column_index() const { return columnIndex_; }; - /** - * @brief Sets the column index to the given value. - * - * @param columnIndex Column index of the cell. - */ - void set_column_index(Index columnIndex) { columnIndex_ = columnIndex; } - - /** - * @brief Assign operator. - */ - Cell_column_index& operator=(Cell_column_index other) { - std::swap(columnIndex_, other.columnIndex_); - return *this; - }; - - private: - Index columnIndex_; /**< Column index. */ -}; - -/** - * @ingroup persistence_matrix - * - * @brief Class managing the value access of a cell. - * - * @tparam Field_element Type of a cell value. - */ -template -class Cell_field_element -{ - public: - /** - * @brief Default constructor. Sets to the element to 0. - */ - Cell_field_element() : element_(0){}; - /** - * @brief Stores the given element. - * - * @param element Value to store. - */ - Cell_field_element(Field_element element) : element_(element){}; - /** - * @brief Copy constructor. - * - * @param cell Cell to copy. - */ - Cell_field_element(const Cell_field_element& cell) : element_(cell.element_){}; - /** - * @brief Move constructor. - * - * @param cell Cell to move. - */ - Cell_field_element(Cell_field_element&& cell) noexcept : element_(std::move(cell.element_)){}; - - /** - * @brief Returns the value stored in the cell. - * - * @return Reference to the value of the cell. - */ - Field_element& get_element() { return element_; }; - /** - * @brief Returns the value stored in the cell. - * - * @return Const reference to the value of the cell. - */ - const Field_element& get_element() const { return element_; }; - /** - * @brief Sets the value. - * - * @param element Value to store. - */ - void set_element(const Field_element& element) { element_ = element; } - - /** - * @brief Assign operator. - */ - Cell_field_element& operator=(Cell_field_element other) { - std::swap(element_, other.element_); - return *this; - }; - - private: - Field_element element_; /**< Value of the cell. */ -}; - -/** - * @class Cell cell_types.h gudhi/Persistence_matrix/columns/cell_types.h - * @ingroup persistence_matrix - * - * @brief %Matrix cell class. Stores by default only the row index it belongs to, but can also store its - * column index when the row access is enabled, as well as its value when they are different from only 0 and 1. - * Zero-valued cells are never made explicit in the matrix. - * - * @tparam Master_matrix An instantiation of @ref Matrix from which all types and options are deduced. - */ -template -class Cell : public Master_matrix::Cell_column_index_option, - public Master_matrix::Cell_field_element_option, - public Master_matrix::Row_hook, - public Master_matrix::Column_hook -{ - private: - using col_opt = typename Master_matrix::Cell_column_index_option; - using field_opt = typename Master_matrix::Cell_field_element_option; - - public: - using Master = Master_matrix; /**< Access to options from outside. */ - using Index = typename Master_matrix::Index; /**< Column index type. */ - using ID_index = typename Master_matrix::ID_index; /**< Row index type. */ - using Field_element = typename Master_matrix::Element; /**< Value type. */ - - /** - * @brief Constructs an cell with all attributes at default values. - */ - Cell(){}; - /** - * @brief Constructs a cell with given row index. Other possible attributes are set at default values. - * - * @param rowIndex @ref rowindex "Row index" of the cell. - */ - Cell(ID_index rowIndex) : col_opt(), field_opt(), rowIndex_(rowIndex){}; - /** - * @brief Constructs a cell with given row and column index. Other possible attributes are set at default values. - * - * @param columnIndex Column index of the cell. - * @param rowIndex @ref rowindex "Row index" of the cell. - */ - Cell(Index columnIndex, ID_index rowIndex) : col_opt(columnIndex), field_opt(), rowIndex_(rowIndex){}; - /** - * @brief Copy constructor. - * - * @param cell Cell to copy. - */ - Cell(const Cell& cell) - : col_opt(static_cast(cell)), - field_opt(static_cast(cell)), - rowIndex_(cell.rowIndex_){}; - /** - * @brief Move constructor. - * - * @param cell Cell to move. - */ - Cell(Cell&& cell) noexcept - : col_opt(std::move(static_cast(cell))), - field_opt(std::move(static_cast(cell))), - rowIndex_(std::exchange(cell.rowIndex_, 0)){}; - - /** - * @brief Returns the row index stored in the cell. - * - * @return @ref rowindex "Row index" of the cell. - */ - ID_index get_row_index() const { return rowIndex_; }; - /** - * @brief Sets the row index stored in the cell. - * - * @param rowIndex @ref rowindex "Row index" of the cell. - */ - void set_row_index(ID_index rowIndex) { rowIndex_ = rowIndex; }; - - /** - * @brief Assign operator. - */ - Cell& operator=(Cell other) { - col_opt::operator=(other); - field_opt::operator=(other); - std::swap(rowIndex_, other.rowIndex_); - return *this; - }; - - /** - * @brief Strictly smaller than comparator. - * - * @param c1 First cell to compare. - * @param c2 Second cell to compare. - * @return true If the row index of the first cell is strictly smaller than the row index of the second cell. - * @return false Otherwise. - */ - friend bool operator<(const Cell& c1, const Cell& c2) { return c1.get_row_index() < c2.get_row_index(); } - /** - * @brief Equality comparator. - * - * @param c1 First cell to compare. - * @param c2 Second cell to compare. - * @return true If the row index of the first cell is equal to the row index of the second cell. - * @return false Otherwise. - */ - friend bool operator==(const Cell& c1, const Cell& c2) { return c1.get_row_index() == c2.get_row_index(); } - - /** - * @brief Converts the cell into a row index. - * - * @return The row index of the cell. - */ - operator ID_index() const { return rowIndex_; } - /** - * @brief Converts the cell into a pair of row index and cell value. - * - * @return A std::pair with first element the row index and second element the value. - */ - operator std::pair() const { - if constexpr (Master_matrix::Option_list::is_z2) { - return {rowIndex_, 1}; - } else { - return {rowIndex_, field_opt::element_}; - } - } - - private: - ID_index rowIndex_; /**< Row index of the cell. */ -}; - -} // namespace persistence_matrix -} // namespace Gudhi - -/** - * @ingroup persistence_matrix - * - * @brief Hash method for @ref Gudhi::persistence_matrix::Cell. - * - * The cells are differentiated by their row indices only. For example, two cells with the same row index - * but different column indices have the same hash value. - * - * @tparam Master_matrix Template parameter of @ref Gudhi::persistence_matrix::Cell. - */ -template -struct std::hash > { - std::size_t operator()(const Gudhi::persistence_matrix::Cell& cell) const { - return std::hash()(cell.get_row_index()); - } -}; - -#endif // PM_MATRIX_CELL_H diff --git a/multipers/gudhi/gudhi/Persistence_matrix/ru_rep_cycles.h b/multipers/gudhi/gudhi/Persistence_matrix/ru_rep_cycles.h index fb2e0c0..e959736 100644 --- a/multipers/gudhi/gudhi/Persistence_matrix/ru_rep_cycles.h +++ b/multipers/gudhi/gudhi/Persistence_matrix/ru_rep_cycles.h @@ -18,9 +18,18 @@ #ifndef PM_RU_REP_CYCLES_H #define PM_RU_REP_CYCLES_H +#include #include //std::move #include //std::sort #include +// #include +#if BOOST_VERSION >= 108100 +#include +#else +#include +#endif + +#include namespace Gudhi { namespace persistence_matrix { @@ -51,9 +60,32 @@ template class RU_representative_cycles { public: - using Index = typename Master_matrix::Index; /**< @ref MatIdx index type. */ - using Bar = typename Master_matrix::Bar; /**< Bar type. */ - using Cycle = typename Master_matrix::Cycle; /**< Cycle type. */ + using Index = typename Master_matrix::Index; /**< @ref MatIdx index type. */ + using ID_index = typename Master_matrix::ID_index; /**< @ref IDIdx index type. */ + using Bar = typename Master_matrix::Bar; /**< Bar type. */ + using Cycle = typename Master_matrix::Cycle; /**< Cycle type. */ + using Cycle_border = std::vector; + using Cycle_borders = std::vector; + struct hashCycle { + size_t operator()(const Cycle_border& cycle) const { + std::hash hasher; + size_t answer = 0; + for (ID_index i : cycle) { + answer ^= hasher(i) + 0x9e3779b9 + (answer << 6) + (answer >> 2); + } + return answer; + } + }; +#if BOOST_VERSION >= 108100 + using Cycle_borders_tmp = boost::unordered_flat_set; +#else + using Cycle_borders_tmp = std::unordered_set; +#endif +#if BOOST_VERSION >= 108100 + using Cycle_unreduced_borders_tmp = boost::unordered_flat_set; +#else + using Cycle_unreduced_borders_tmp = std::unordered_set; +#endif /** * @brief Default constructor. @@ -77,22 +109,30 @@ class RU_representative_cycles */ void update_representative_cycles(); + // /** + // * @brief Returns the current representative cycles. If the matrix is modified later after the first call, + // * @ref update_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(); + // /** + // * @brief Returns the representative cycle corresponding to the given bar. + // * If the matrix is modified later after the first call, + // * @ref update_representative_cycles has to be called to update the returned cycles. + // * + // * @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); + /** * @brief Returns the current representative cycles. If the matrix is modified later after the first call, * @ref update_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(); - /** - * @brief Returns the representative cycle corresponding to the given bar. - * If the matrix is modified later after the first call, - * @ref update_representative_cycles has to be called to update the returned cycles. - * - * @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); + std::vector get_representative_cycles_as_borders(bool detailed = false); /** * @brief Assign operator. @@ -103,91 +143,137 @@ class RU_representative_cycles */ friend void swap(RU_representative_cycles& base1, RU_representative_cycles& base2) { base1.representativeCycles_.swap(base2.representativeCycles_); - base1.birthToCycle_.swap(base2.birthToCycle_); + base1.u_transposed_.swap(base2.u_transposed_); + // base1.birthToCycle_.swap(base2.birthToCycle_); } private: using Master_RU_matrix = typename Master_matrix::Master_RU_matrix; - std::vector representativeCycles_; /**< Cycle container. */ - std::vector birthToCycle_; /**< Map from birth index to cycle index. */ + std::vector representativeCycles_; /**< Cycle container. */ + // std::vector birthToCycle_; /**< Map from birth index to cycle index. */ + std::vector u_transposed_; + + void _get_initial_borders(Index idx, Cycle_borders_tmp& borders); + Cycle_border _get_border(Index uIndex); + Cycle_border _get_dependent_border(Index uIndex); constexpr Master_RU_matrix* _matrix() { return static_cast(this); } constexpr const Master_RU_matrix* _matrix() const { return static_cast(this); } }; template -inline RU_representative_cycles::RU_representative_cycles() +inline RU_representative_cycles::RU_representative_cycles() {} template inline RU_representative_cycles::RU_representative_cycles( const RU_representative_cycles& matrixToCopy) - : representativeCycles_(matrixToCopy.representativeCycles_), birthToCycle_(matrixToCopy.birthToCycle_) + : representativeCycles_(matrixToCopy.representativeCycles_), + // birthToCycle_(matrixToCopy.birthToCycle_), + u_transposed_(matrixToCopy.u_transposed_) {} template inline RU_representative_cycles::RU_representative_cycles( RU_representative_cycles&& other) noexcept - : representativeCycles_(std::move(other.representativeCycles_)), birthToCycle_(std::move(other.birthToCycle_)) + : representativeCycles_(std::move(other.representativeCycles_)), + // birthToCycle_(std::move(other.birthToCycle_)), + u_transposed_(std::move(other.u_transposed_)) {} +//TODO: u_transposed_ as a vector of Index and another vector giving position + length of each cycle? template inline void RU_representative_cycles::update_representative_cycles() { - if constexpr (Master_matrix::Option_list::is_z2) { - birthToCycle_.clear(); - birthToCycle_.resize(_matrix()->reducedMatrixR_.get_number_of_columns(), -1); - Index c = 0; - for (Index i = 0; i < _matrix()->reducedMatrixR_.get_number_of_columns(); i++) { - if (_matrix()->reducedMatrixR_.is_zero_column(i)) { - birthToCycle_[i] = c; - ++c; - } + //WARNING: this was only thought for the multipers interface, it is not definitive and does not + // cover all cases. + static_assert(Master_matrix::Option_list::has_column_pairings && Master_matrix::Option_list::is_z2, + "Needs an ID to Pos map."); + + auto rsize = _matrix()->reducedMatrixR_.get_number_of_columns(); + representativeCycles_.clear(); + representativeCycles_.reserve(rsize); + + auto usize = _matrix()->mirrorMatrixU_.get_number_of_columns(); + u_transposed_.clear(); + u_transposed_.resize(usize); + + for (Index i = 0; i < usize; i++) { + if (_matrix()->reducedMatrixR_.is_zero_column(i)) { + representativeCycles_.push_back(i); } - representativeCycles_.clear(); - representativeCycles_.resize(c); - for (Index i = 0; i < _matrix()->mirrorMatrixU_.get_number_of_columns(); i++) { - for (const auto& entry : _matrix()->mirrorMatrixU_.get_column(i)) { - auto idx = birthToCycle_[entry.get_row_index()]; - if (idx != static_cast(-1)) { - representativeCycles_[idx].push_back(i); + if constexpr (Master_matrix::Option_list::column_type == Column_types::HEAP || + Master_matrix::Option_list::column_type == Column_types::VECTOR) { + // TODO: have a better way to do this. For now, one cannot use the column iterators for that. + unsigned int j = 0; + for (const auto& entry : _matrix()->mirrorMatrixU_.get_column(i).get_content()) { + if (entry != 0){ + u_transposed_[j].push_back(i); } + ++j; } - } - } else { - birthToCycle_.clear(); - birthToCycle_.resize(_matrix()->reducedMatrixR_.get_number_of_columns(), -1); - for (Index i = 0; i < _matrix()->reducedMatrixR_.get_number_of_columns(); i++) { - if (_matrix()->reducedMatrixR_.is_zero_column(i)) { - representativeCycles_.push_back(Cycle()); - for (const auto& entry : _matrix()->mirrorMatrixU_.get_column(i)) { - representativeCycles_.back().push_back(entry.get_row_index()); - } - if constexpr (std::is_same_v || - std::is_same_v) - std::sort(representativeCycles_.back().begin(), representativeCycles_.back().end()); - birthToCycle_[i] = representativeCycles_.size() - 1; + } else { + for (const auto& entry : _matrix()->mirrorMatrixU_.get_column(i)) { + u_transposed_[entry.get_row_index()].push_back(i); } } } } -template -inline const std::vector::Cycle>& -RU_representative_cycles::get_representative_cycles() -{ - if (representativeCycles_.empty()) update_representative_cycles(); - return representativeCycles_; -} +// template +// inline const std::vector::Cycle>& +// RU_representative_cycles::get_representative_cycles() +// { +// if (representativeCycles_.empty()) update_representative_cycles(); +// return representativeCycles_; +// } + +// template +// inline const typename RU_representative_cycles::Cycle& +// RU_representative_cycles::get_representative_cycle(const Bar& bar) +// { +// if (representativeCycles_.empty()) update_representative_cycles(); +// return representativeCycles_[birthToCycle_[bar.birth]]; +// } +//TODO: not store cycle borders in vectors of vectors of vectors template -inline const typename RU_representative_cycles::Cycle& -RU_representative_cycles::get_representative_cycle(const Bar& bar) +inline std::vector::Cycle_borders> +RU_representative_cycles::get_representative_cycles_as_borders(bool detailed) { + static_assert(Master_matrix::Option_list::is_z2, "Only available for Z2 coefficients for now."); + if (representativeCycles_.empty()) update_representative_cycles(); - return representativeCycles_[birthToCycle_[bar.birth]]; + + std::vector cycles(representativeCycles_.size()); + unsigned int i = 0; + for (const auto& cycleIndex : representativeCycles_){ + const auto& cycle = u_transposed_[cycleIndex]; + assert(_matrix()->reducedMatrixR_.is_zero_column(cycle.back())); + if (cycle.size() != 1){ //cycle.size == 1 -> border is empty + if (detailed) { + Cycle_borders_tmp cbt; + for (unsigned int j = 0; j < cycle.size() - 1; ++j){ + assert(!_matrix()->reducedMatrixR_.is_zero_column(cycle[j])); + _get_initial_borders(cycle[j], cbt); + } + cycles[i] = Cycle_borders(cbt.begin(), cbt.end()); + cycles[i].push_back(_get_dependent_border(cycle.back())); + } else { + cycles[i].resize(cycle.size()); + for (unsigned int j = 0; j < cycle.size() - 1; ++j){ + cycles[i][j] = _get_border(cycle[j]); + } + cycles[i].back() = _get_dependent_border(cycle.back()); + } + } else { + cycles[i].push_back({}); + } + + ++i; + } + return cycles; } template @@ -195,10 +281,95 @@ inline RU_representative_cycles& RU_representative_cycles other) { representativeCycles_.swap(other.representativeCycles_); - birthToCycle_.swap(other.birthToCycle_); + u_transposed_.swap(other.u_transposed_); + // birthToCycle_.swap(other.birthToCycle_); return *this; } +//TODO: try avoid storing the vectors in the unordered set, when it is not reduced as it adds run time for hash +//computation and possibly comparaison of vectors. Treat them separately. See if it goes faster, or if it was marginal +template +inline void RU_representative_cycles::_get_initial_borders(Index idx, + Cycle_borders_tmp& borders) { + const Cycle& cycle = u_transposed_[idx]; + assert(cycle.back() == idx); + + auto add_to = [](const std::vector& b, Cycle_borders_tmp& borders){ + auto p = borders.insert(b); + if (!p.second) borders.erase(p.first); + }; + + for (unsigned int j = 0; j < cycle.size() - 1; ++j){ + _get_initial_borders(cycle[j], borders); + } + + add_to(_get_dependent_border(idx), borders); +} + +template +inline RU_representative_cycles::Cycle_border +RU_representative_cycles::_get_border(Index uIndex) { + const auto& col = _matrix()->reducedMatrixR_.get_column(uIndex); + Cycle_border res; + + unsigned int j = 0; + if constexpr (Master_matrix::Option_list::column_type == Column_types::HEAP || + Master_matrix::Option_list::column_type == Column_types::VECTOR) { + res.reserve(col.size()); + // TODO: have a better way to do this. For now, one cannot use the column iterators for that. + for (auto i : col.get_content()) { + if (i != 0) res.push_back(j); + ++j; + } + } else { + res.resize(col.size()); + for (const auto& entry : col) { + res[j] = entry.get_row_index(); + ++j; + } + } + + return res; +} + +template +inline RU_representative_cycles::Cycle_border +RU_representative_cycles::_get_dependent_border(Index uIndex) { + if (u_transposed_[uIndex].size() == 1){ + return _get_border(uIndex); + } + + auto add = [](const Master_matrix::Column& col, Cycle_unreduced_borders_tmp& b){ + if constexpr (Master_matrix::Option_list::column_type == Column_types::HEAP || + Master_matrix::Option_list::column_type == Column_types::VECTOR) { + // TODO: have a better way to do this. For now, one cannot use the column iterators for that. + unsigned int j = 0; + for (auto k : col.get_content()) { + if (k != 0) { + auto p = b.insert(j); + if (!p.second) b.erase(p.first); + } + ++j; + } + } else { + for (const auto& entry : col) { + auto p = b.insert(entry.get_row_index()); + if (!p.second) b.erase(p.first); + } + } + }; + + Cycle_unreduced_borders_tmp b; + for (Index i : u_transposed_[uIndex]){ + add(_matrix()->reducedMatrixR_.get_column(i), b); + } + + Cycle_border res(b.begin(), b.end()); + std::sort(res.begin(), res.end()); + + return res; +} + } // namespace persistence_matrix } // namespace Gudhi diff --git a/multipers/gudhi/mma_interface_matrix.h b/multipers/gudhi/mma_interface_matrix.h index c1c2cef..12acefe 100644 --- a/multipers/gudhi/mma_interface_matrix.h +++ b/multipers/gudhi/mma_interface_matrix.h @@ -23,7 +23,6 @@ #include #include -#include #include #include @@ -74,6 +73,7 @@ class Persistence_backend_matrix { // using index = typename matrix_type::index; // using id_index = typename matrix_type::id_index; using pos_index = typename matrix_type::Pos_index; + using Index = typename matrix_type::Index; using dimension_type = typename matrix_type::Dimension; class Barcode_iterator @@ -263,28 +263,37 @@ class Persistence_backend_matrix { return stream; } - inline std::vector get_representative_cycles(bool update) { - // Only used when vineyard, so shrinked permutation i.e. + inline std::vector>> get_representative_cycles(bool update, + bool detailed = false) { + // Only used when vineyard, so shrunk permutation i.e. // without the -1, is permutation as we keep inf values (they can become // finite) cf barcode perm which is copied to remove the -1 + std::vector permutation2; + permutation2.reserve(permutation_->size()); + for (std::size_t i : *permutation_) { + if (i >= matrix_.get_number_of_columns()) { + continue; + } + permutation2.push_back(i); + } const bool verbose = false; if (update) [[likely]] matrix_.update_representative_cycles(); - auto current_cycles = matrix_.get_representative_cycles(); - for (auto &truc : current_cycles) { - if constexpr (verbose) - std::cout << "Cycle (matrix_ order): "; - for (auto &machin : truc) { - if constexpr (verbose) { - std::cout << " matrix order: " << machin; - std::cout << ", structure order : " << permutation_[machin] - << std::endl; + std::vector>> current_cycles = + matrix_.get_representative_cycles_as_borders(detailed); + unsigned int i = 0; + for (auto &cycle : current_cycles) { + if constexpr (verbose) std::cout << "Cycle (matrix_ order): "; + for (auto &border : cycle) { + for (auto &b : border) { + b = permutation2[b]; } - machin = permutation_->operator[](machin); } + ++i; } return current_cycles; } + inline void _update_permutation_ptr(std::vector &perm) { permutation_ = &perm; } diff --git a/multipers/gudhi/truc.h b/multipers/gudhi/truc.h index ee713db..5df267a 100644 --- a/multipers/gudhi/truc.h +++ b/multipers/gudhi/truc.h @@ -647,35 +647,19 @@ class Truc { } // dim, num_cycle_of_dim, num_faces_in_cycle, vertices_in_face - inline std::vector>>> - get_representative_cycles(bool update = true, - const std::set &dims = {}) { + inline std::vector>>> get_representative_cycles( + bool update = true, bool detailed = false) { // iterable iterable simplex key - auto cycles_key = persistence.get_representative_cycles(update); + auto cycles_key = persistence.get_representative_cycles(update, detailed); auto num_cycles = cycles_key.size(); - std::vector>>> out( - structure.max_dimension() + 1); - for (auto &cycles_of_dim : out) - cycles_of_dim.reserve(num_cycles); - for (auto &cycle : cycles_key) { - if (structure.dimension(cycle[0]) == 0 && - (dims.size() == 0 || dims.contains(0))) { - out[0].push_back({{static_cast(cycle[0])}}); - continue; + std::vector>>> out(structure.max_dimension() + 1); + for (auto &cycles_of_dim : out) cycles_of_dim.reserve(num_cycles); + for (const auto &cycle : cycles_key) { + int cycle_dim = 0; // for more generality, should be minimal dimension instead + if (!cycle[0].empty()) { // if empty, cycle has no border -> assumes dimension 0 even if it could be min dim + cycle_dim = structure.dimension(cycle[0][0]) + 1; // all faces have the same dim } - int cycle_dim = - structure.dimension(cycle[0]); // all faces have the same dim - // if (!(dims.size() == 0 || dims.contains(cycle_dim))) - // continue; - std::vector> cycle_faces(cycle.size()); - for (std::size_t i = 0; i < cycle_faces.size(); ++i) { - const auto &B = structure[cycle[i]]; - cycle_faces[i].resize(B.size()); - for (std::size_t j = 0; j < B.size(); ++j) { - cycle_faces[i][j] = static_cast(B[j]); - } - } - out[cycle_dim].push_back(cycle_faces); + out[cycle_dim].push_back(cycle); } return out; } diff --git a/multipers/slicer.pxd.tp b/multipers/slicer.pxd.tp index 6534562..0ca19eb 100644 --- a/multipers/slicer.pxd.tp +++ b/multipers/slicer.pxd.tp @@ -76,7 +76,7 @@ cdef extern from "Persistence_slices_interface.h": {{endif}} {{if IS_VINE}} void vineyard_update() nogil - vector[vector[vector[vector[unsigned int]]]] get_representative_cycles(bool) nogil + vector[vector[vector[vector[unsigned int]]]] get_representative_cycles(bool, bool) nogil vector[size_t] get_current_order() nogil {{endif}} diff --git a/multipers/slicer.pyx.tp b/multipers/slicer.pyx.tp index 05531c6..2ad5f7c 100644 --- a/multipers/slicer.pyx.tp +++ b/multipers/slicer.pyx.tp @@ -438,12 +438,12 @@ cdef class {{PYTHON_TYPE}}: self.push_to_line(basepoint,direction) self.truc.vineyard_update() return self - def get_representative_cycles(self, bool update=True): + def get_representative_cycles(self, bool update=True, bool detailed=False): """ Returns the representative cycles of the current barcode. Recomputes the generators if update=True """ - return self.truc.get_representative_cycles(update) + return self.truc.get_representative_cycles(update, detailed) def get_permutation(self): """ Returns the current generator permutation (w.r.t. vineyard). diff --git a/multipers/tests/test_slicer.py b/multipers/tests/test_slicer.py index b34b6ab..2944a75 100644 --- a/multipers/tests/test_slicer.py +++ b/multipers/tests/test_slicer.py @@ -96,7 +96,7 @@ def test_representative_cycles(): cycles = slicer.get_representative_cycles() assert len(cycles) == 3, f"There should be 3 dimensions here, found {len(cycles)}" assert ( - np.asarray(cycles[0]).size == 7 + np.asarray(cycles[0]).shape[0] == 7 ), f"Invalid number of 0-cycles, got {np.asarray(cycles[0]).size}" for c in cycles[1]: assert (