From 353929a1ff36d159ccf797a2adeed71761845165 Mon Sep 17 00:00:00 2001 From: hschreiber Date: Thu, 4 Apr 2024 19:23:43 +0200 Subject: [PATCH] doc --- .../concept/FieldOperators.h | 3 +- .../concept/PersistenceMatrixColumn.h | 452 ++++ .../concept/PersistenceMatrixOptions.h | 3 +- .../columns/cell_constructors.h | 149 +- .../Persistence_matrix/columns/cell_types.h | 457 ++-- .../columns/chain_column_extra_properties.h | 163 +- .../columns/column_dimension_holder.h | 109 +- .../Persistence_matrix/columns/heap_column.h | 1903 +++++++++-------- .../columns/intrusive_list_column.h | 1822 ++++++++-------- .../columns/intrusive_set_column.h | 1778 +++++++-------- src/Persistence_matrix/include/gudhi/matrix.h | 5 +- 11 files changed, 3890 insertions(+), 2954 deletions(-) create mode 100644 src/Persistence_matrix/concept/PersistenceMatrixColumn.h diff --git a/src/Persistence_matrix/concept/FieldOperators.h b/src/Persistence_matrix/concept/FieldOperators.h index 47ed2a0d99..1f5173e51d 100644 --- a/src/Persistence_matrix/concept/FieldOperators.h +++ b/src/Persistence_matrix/concept/FieldOperators.h @@ -22,7 +22,8 @@ namespace persistence_matrix { * Implementations of this concept are @ref Zp_field_operators, @ref Z2_field_operators, * @ref Multi_field_operators and @ref Multi_field_small_operators. */ -struct FieldOperators { +class FieldOperators +{ public: using element_type = unspecified; /**< Type for the elements in the field. */ using characteristic_type = unspecified; /**< Type for the field characteristic. */ diff --git a/src/Persistence_matrix/concept/PersistenceMatrixColumn.h b/src/Persistence_matrix/concept/PersistenceMatrixColumn.h new file mode 100644 index 0000000000..dd97f01bab --- /dev/null +++ b/src/Persistence_matrix/concept/PersistenceMatrixColumn.h @@ -0,0 +1,452 @@ +/* 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 PersistenceMatrixColumn.h + * @author Hannah Schreiber + * @brief Contains the concept for the matrix columns. + */ + +namespace Gudhi { +namespace persistence_matrix { + +/** @brief Concept of the column classes used by the @ref Matrix class. + * + * Implementations of this concept are @ref Heap_column, @ref List_column, @ref Vector_column, @ref Naive_vector_column + * @ref Set_column, @ref Unordered_set_column, @ref PersistenceMatrixColumn and @ref Intrusive_set_column. + */ +class PersistenceMatrixColumn : + /** + * @brief If PersistenceMatrixOptions::has_row_access is true, then @ref Row_access. + * Otherwise @ref Dummy_row_access. Can eventually be removed if the structure of the column does not allow + * row access (as for @ref Heap_column), but then it needs to be notified in the documentation of @ref Column_types + * and as static_assert in @ref Matrix::_assert_options. + */ + public Row_access_option, + /** + * @brief If PersistenceMatrixOptions::has_column_pairings || + PersistenceMatrixOptions::has_vine_update || + PersistenceMatrixOptions::can_retrieve_representative_cycles is true, then @ref + Column_dimension_holder. Otherwise @ref Dummy_dimension_holder. + * + */ + public Column_dimension_option, + /** + * @brief If (PersistenceMatrixOptions::has_column_pairings || + PersistenceMatrixOptions::has_vine_update || + PersistenceMatrixOptions::can_retrieve_representative_cycles) && + !PersistenceMatrixOptions::is_of_boundary_type is true, then @ref Chain_column_extra_properties. + Otherwise @ref Dummy_chain_properties. + */ + public Chain_column_option +{ + public: + using Master = unspecified; /**< Master matrix. */ + using Cell_constructor = unspecified; /**< @ref Cell factory. */ + using Field_operators = unspecified; /**< Follows the @ref FieldOperators concept. */ + using Field_element_type = unspecified; /**< Type of a field element. */ + using index = unspecified; /**< Type of MatIdx index. */ + using id_index = unspecified; /**< Type of IDIdx index. */ + using dimension_type = unspecified; /**< Type for dimension value. */ + using Cell = unspecified; /**< @ref Cell. */ + using Column_type = unspecified; /**< Type of cell container. */ + using iterator = unspecified; /**< Column iterator type. */ + using const_iterator = unspecified; /**< Column const_iterator type. */ + using reverse_iterator = unspecified; /**< Column reverse_iterator type. */ + using const_reverse_iterator = unspecified; /**< Column const_reverse_iterator type. */ + + /** + * @brief Constructs an empty column. If @p cellConstructor is not specified or is set to nullptr, the column + * can only be used as a dummy, i.e., no modifying method should be used or there will be a segmentation fault. + * Same goes for @p operators if @ref PersistenceMatrixOptions::is_z2 is false. + * + * @param operators Pointer to the field operators. + * @param cellConstructor Pointer to the cell factory. + */ + PersistenceMatrixColumn(Field_operators* operators = nullptr, Cell_constructor* cellConstructor = nullptr); + /** + * @brief Constructs a column from the given range of @ref Matrix::cell_rep_type. If the dimension is stored, + * the face is assumed to be simplicial and its dimension to be @ref nonZeroRowIndices length - 1 or 0. + * Otherwise, the dimension should be specified with another constructor. + * + * @tparam Container_type Range of @ref Matrix::cell_rep_type. Assumed to have a begin(), end() and size() method. + * @param nonZeroRowIndices Range of @ref Matrix::cell_rep_type representing all rows with non zero values. + * @param operators Pointer to the field operators. + * @param cellConstructor Pointer to the cell factory. + */ + template + PersistenceMatrixColumn(const Container_type& nonZeroRowIndices, + Field_operators* operators, + Cell_constructor* cellConstructor); + /** + * @brief Constructs a column from the given range of @ref Matrix::cell_rep_type such that the rows can be accessed. + * Each new cell in the column is also inserted in a row using @ref Row_access::insert_cell. + * If the dimension is stored, the face is assumed to be simplicial and its dimension to be + * @ref nonZeroRowIndices length - 1 or 0. Otherwise, the dimension should be specified with another constructor. + * + * @tparam Container_type Range of @ref Matrix::cell_rep_type. Assumed to have a begin(), end() and size() method. + * @tparam Row_container_type Either std::map if @ref PersistenceMatrixOptions::has_removable_rows is true or + * std::vector. + * @param columnIndex MatIdx column index that should be specified to the cells. + * @param nonZeroRowIndices Range of @ref Matrix::cell_rep_type representing all rows with non zero values. + * @param rowContainer Pointer to the row container that will be forwarded to @ref Row_access at construction. + * @param operators Pointer to the field operators. + * @param cellConstructor Pointer to the cell factory. + */ + template + PersistenceMatrixColumn(index columnIndex, + const Container_type& nonZeroRowIndices, + Row_container_type* rowContainer, + Field_operators* operators, + Cell_constructor* cellConstructor); + /** + * @brief Constructs a column from the given range of @ref Matrix::cell_rep_type and stores the given dimension + * if @ref Column_dimension_option is not a dummy. + * + * @tparam Container_type Range of @ref Matrix::cell_rep_type. Assumed to have a begin(), end() and size() method. + * @param nonZeroChainRowIndices Range of @ref Matrix::cell_rep_type representing all rows with non zero values. + * @param dimension Dimension of the column. Is ignored if the dimension is not stored. + * @param operators Pointer to the field operators. + * @param cellConstructor Pointer to the cell factory. + */ + template + PersistenceMatrixColumn(const Container_type& nonZeroChainRowIndices, + dimension_type dimension, + Field_operators* operators, + Cell_constructor* cellConstructor); + /** + * @brief Constructs a column from the given range of @ref Matrix::cell_rep_type such that the rows can be accessed. + * Each new cell in the column is also inserted in a row using @ref Row_access::insert_cell. + * Stores the given dimension if @ref Column_dimension_option is not a dummy. + * + * @tparam Container_type Range of @ref Matrix::cell_rep_type. Assumed to have a begin(), end() and size() method. + * @tparam Row_container_type Either std::map if @ref PersistenceMatrixOptions::has_removable_rows is true or + * std::vector. + * @param columnIndex MatIdx column index that should be specified to the cells. + * @param nonZeroRowIndices Range of @ref Matrix::cell_rep_type representing all rows with non zero values. + * @param dimension Dimension of the column. Is ignored if the dimension is not stored. + * @param rowContainer Pointer to the row container that will be forwarded to @ref Row_access at construction. + * @param operators Pointer to the field operators. + * @param cellConstructor Pointer to the cell factory. + */ + template + PersistenceMatrixColumn(index columnIndex, + const Container_type& nonZeroChainRowIndices, + dimension_type dimension, + Row_container_type* rowContainer, + Field_operators* operators, + Cell_constructor* cellConstructor); + /** + * @brief Copy constructor. If @p operators or @p cellConstructor is not a null pointer, its value is kept + * instead of the one in the copied column. + * + * @param column Column to copy. + * @param operators Pointer to the field operators. + * If null pointer, the pointer in @p column is choosen instead. + * @param cellConstructor Pointer to the cell factory. + * If null pointer, the pointer in @p column is choosen instead. + */ + PersistenceMatrixColumn(const PersistenceMatrixColumn& column, + Field_operators* operators = nullptr, + Cell_constructor* cellConstructor = nullptr); + /** + * @brief Copy constructor with row access. + * If @p operators or @p cellConstructor is not a null pointer, its value is kept + * instead of the one in the copied column. + * + * @tparam Row_container_type Either std::map if @ref PersistenceMatrixOptions::has_removable_rows is true or + * std::vector. + * @param column Column to copy. + * @param columnIndex MatIdx column index of the new column once copied. + * @param rowContainer Pointer to the row container that will be forwarded to @ref Row_access. + * @param operators Pointer to the field operators. + * If null pointer, the pointer in @p column is choosen instead. + * @param cellConstructor Pointer to the cell factory. + * If null pointer, the pointer in @p column is choosen instead. + */ + template + PersistenceMatrixColumn(const PersistenceMatrixColumn& column, + index columnIndex, + Row_container_type* rowContainer, + Field_operators* operators = nullptr, + Cell_constructor* cellConstructor = nullptr); + /** + * @brief Move constructor. + * + * @param column Column to move. + */ + PersistenceMatrixColumn(PersistenceMatrixColumn&& column) noexcept; + /** + * @brief Destructor. + */ + ~PersistenceMatrixColumn(); + + /** + * @brief Returns the values of the column, zero values included. + * + * @param columnLength Number of rows to be returned. If -1, the number of rows is fixed at the biggest + * row index with non zero value. Default value: -1. + * @return Vector of @ref Field_element_type. At element \f$ i \f$ of the vector will be stored the value + * at row \f$ i \f$ of the column. + */ + std::vector get_content(int columnLength = -1) const; + /** + * @brief Indicates if the cell at given row index has value zero. + * + * @param rowIndex Row index to look at. + * @return true If the cell has value zero. + * @return false Otherwise. + */ + bool is_non_zero(id_index rowIndex) const; + /** + * @brief Indicates if the column is empty or has only zero values. + * + * @return true If the column is empty or only has zero values. + * @return false Otherwise. + */ + bool is_empty(); + /** + * @brief Returns the size of the underlying container. + * + * @warning Depending of the column type, the container does not have to contain only the non-zero cells. + * Even if for most of the types, the size of the container will correspond to the number of non-zero values + * in the column, it is not always the case. See description of the actual Column class for more details. + * + * @return Size of the underlying container. + */ + std::size_t size() const; + + /** + * @brief Reorders the column with the given map of row indices. Also changes the column index stored in the + * cells if row access is enabled and @p columnIndex is not -1. + * + * Only useful for base and boundary matrices using lazy swaps. + * + * @tparam Map_type Map with an `at` method. + * @param valueMap Map such that `valueMap.at(i)` indicates the new row index of the cell + * at current row index `i`. + * @param columnIndex New MatIdx column index of the column. If -1, the index does not change. Ignored if + * the row access is not enabled. Default value: -1. + */ + template + void reorder(const Map_type& valueMap, [[maybe_unused]] index columnIndex = -1); + /** + * @brief Zeros/empties the column. + * + * Only useful for base and boundary matrices. Used in `zero_column` and in the reduction algorithm + * for the persistence barcode. + */ + void clear(); + /** + * @brief Zeros the cell at given row index. + * + * Only useful for base and boundary matrices. Used in `zero_cell` and during vine swaps. + * + * @param rowIndex Row index of the cell to zero. + */ + void clear(id_index rowIndex); + + /** + * @brief Returns the row index of the pivot. If the column does not have a pivot, returns -1. + * + * Only useful for boundary and chain matrices. + * + * @return Row index of the pivot or -1. + */ + id_index get_pivot(); + /** + * @brief Returns the value of the pivot. If the column does not have a pivot, returns 0. + * + * Has to have value 1 if \f$ Z_2 \f$ coefficients are used. + * + * Only useful for boundary and chain matrices. + * + * @return The value of the pivot or 0. + */ + Field_element_type get_pivot_value(); + + /** + * @brief Returns a begin @ref Cell iterator to iterate over all cells contained in the underlying container. + * + * @warning The iterators really just iterate over the underlying container. Depending of the column type, + * neither the content nor the order is garanteed. See description of the actual Column class for more details. + * + * @return @ref Cell iterator. + */ + iterator begin() noexcept; + /** + * @brief Returns a begin @ref Cell const iterator to iterate over all cells contained in the underlying container. + * + * @warning The iterators really just iterate over the underlying container. Depending of the column type, + * neither the content nor the order is garanteed. See description of the actual Column class for more details. + * + * @return @ref Cell const iterator. + */ + const_iterator begin() const noexcept; + /** + * @brief Returns a end @ref Cell iterator, iterating over all cells contained in the underlying container. + * + * @warning The iterators really just iterate over the underlying container. Depending of the column type, + * neither the content nor the order is garanteed. See description of the actual Column class for more details. + * + * @return @ref Cell iterator. + */ + iterator end() noexcept; + /** + * @brief Returns a end @ref Cell const iterator, iterating over all cells contained in the underlying container. + * + * @warning The iterators really just iterate over the underlying container. Depending of the column type, + * neither the content nor the order is garanteed. See description of the actual Column class for more details. + * + * @return @ref Cell const iterator. + */ + const_iterator end() const noexcept; + /** + * @brief Returns a begin @ref Cell reverse iterator to iterate over all cells contained in the underlying container. + * + * @warning The iterators really just iterate over the underlying container. Depending of the column type, + * neither the content nor the order is garanteed. See description of the actual Column class for more details. + * + * @return @ref Cell reverse iterator. + */ + reverse_iterator rbegin() noexcept; + /** + * @brief Returns a begin @ref Cell const reverse iterator to iterate over all cells contained in the underlying + * container. + * + * @warning The iterators really just iterate over the underlying container. Depending of the column type, + * neither the content nor the order is garanteed. See description of the actual Column class for more details. + * + * @return @ref Cell const reverse iterator. + */ + const_reverse_iterator rbegin() const noexcept; + /** + * @brief Returns a end @ref Cell reverse iterator, iterating over all cells contained in the underlying container. + * + * @warning The iterators really just iterate over the underlying container. Depending of the column type, + * neither the content nor the order is garanteed. See description of the actual Column class for more details. + * + * @return @ref Cell reverse iterator. + */ + reverse_iterator rend() noexcept; + /** + * @brief Returns a end @ref Cell const reverse iterator, iterating over all cells contained in the underlying + * container. + * + * @warning The iterators really just iterate over the underlying container. Depending of the column type, + * neither the content nor the order is garanteed. See description of the actual Column class for more details. + * + * @return @ref Cell const reverse iterator. + */ + const_reverse_iterator rend() const noexcept; + + /** + * @brief Adds the given @ref Cell range onto the column. + * + * @tparam Cell_range @ref Cell range with `begin` and `end` method. + * @param column @ref Cell range. Every cell has to return the right value when using `get_row_index` and, + * if @ref PersistenceMatrixOptions::is_z2 is false, also when using `get_element`. + * @return Reference to this column. + */ + template + PersistenceMatrixColumn& operator+=(const Cell_range& column); + /** + * @brief Adds the given column onto this column. + * + * @param column Column to add. + * @return Reference to this column. + */ + PersistenceMatrixColumn& operator+=(PersistenceMatrixColumn& column); + + /** + * @brief Multiplies all values in the column with the given value. + * + * @param val Value to multiply. + * @return Reference to this column. + */ + PersistenceMatrixColumn& operator*=(const Field_element_type& val); + + /** + * @brief @p this = @p val * @p this + @p column + * + * @tparam Cell_range @ref Cell range with `begin` and `end` method. + * @param val Value to multiply. + * @param column @ref Cell range. Every cell has to return the right value when using `get_row_index` and, + * if @ref PersistenceMatrixOptions::is_z2 is false, also when using `get_element`. + * @return Reference to this column. + */ + template + PersistenceMatrixColumn& multiply_and_add(const Field_element_type& val, const Cell_range& column); + /** + * @brief @p this = @p val * @p this + @p column + * + * @param val Value to multiply. + * @param column Column to add. + * @return Reference to this column. + */ + PersistenceMatrixColumn& multiply_and_add(const Field_element_type& val, PersistenceMatrixColumn& column); + /** + * @brief @p this = @p this + @p column * @p val + * + * @tparam Cell_range @ref Cell range with `begin` and `end` method. + * @param column @ref Cell range. Every cell has to return the right value when using `get_row_index` and, + * if @ref PersistenceMatrixOptions::is_z2 is false, also when using `get_element`. + * @param val Value to multiply. + * @return Reference to this column. + */ + template + PersistenceMatrixColumn& multiply_and_add(const Cell_range& column, const Field_element_type& val); + /** + * @brief @p this = @p this + @p column * @p val + * + * @param column Column to add. + * @param val Value to multiply. + * @return Reference to this column. + */ + PersistenceMatrixColumn& multiply_and_add(PersistenceMatrixColumn& column, const Field_element_type& val); + + /** + * @brief Equality comparator. Equal in the sense that what is "supposed" to be contained in the columns is equal, + * not what is actually stored in the underlying container. For exemple, the underlying container of + * @ref Vector_column can contain cells which were erased explicitely by @ref clear(index). Those cells should not + * be taken into account while comparing. + * + * @param c1 First column to compare. + * @param c2 Second column to compare. + * @return true If both column are equal. + * @return false Otherwise. + */ + friend bool operator==(const PersistenceMatrixColumn& c1, const PersistenceMatrixColumn& c2); + /** + * @brief "Strictly smaller than" comparator. Usually a lexicographical order, but what matters is that the + * order is total. The order should apply on what is "supposed" to be contained in the columns, + * not what is actually stored in the underlying container. For exemple, the underlying container of + * @ref Vector_column can contain cells which were erased explicitely by @ref clear(index). Those cells should not + * be taken into account while comparing. + * + * @param c1 First column to compare. + * @param c2 Second column to compare. + * @return true If the first column is strictly smaller than the second one. + * @return false Otherwise. + */ + friend bool operator<(const PersistenceMatrixColumn& c1, const PersistenceMatrixColumn& c2); + + /** + * @brief Assign operator. Should be disabled when row access is enabled. + */ + PersistenceMatrixColumn& operator=(const PersistenceMatrixColumn& other); + /** + * @brief Swap operator. + */ + friend void swap(PersistenceMatrixColumn& col1, PersistenceMatrixColumn& col2); +}; + +} // namespace persistence_matrix +} // namespace Gudhi diff --git a/src/Persistence_matrix/concept/PersistenceMatrixOptions.h b/src/Persistence_matrix/concept/PersistenceMatrixOptions.h index ed2ac56536..97e540e52f 100644 --- a/src/Persistence_matrix/concept/PersistenceMatrixOptions.h +++ b/src/Persistence_matrix/concept/PersistenceMatrixOptions.h @@ -23,7 +23,8 @@ namespace persistence_matrix { * If you want to provide your own, it is recommended that you derive from it and override some parts instead of * writing a class from scratch. */ -struct PersistenceMatrixOptions { +struct PersistenceMatrixOptions +{ /** * @brief Field operators. Has to follow the @ref [TODO: concept] concept. * The type will not be used if @ref is_z2 is set to true, so it can be set to anything. diff --git a/src/Persistence_matrix/include/gudhi/Persistence_matrix/columns/cell_constructors.h b/src/Persistence_matrix/include/gudhi/Persistence_matrix/columns/cell_constructors.h index 532cf51d78..160189f8ee 100644 --- a/src/Persistence_matrix/include/gudhi/Persistence_matrix/columns/cell_constructors.h +++ b/src/Persistence_matrix/include/gudhi/Persistence_matrix/columns/cell_constructors.h @@ -8,58 +8,125 @@ * - YYYY/MM Author: Description of the modification */ +/** + * @file cell_constructors.h + * @author Hannah Schreiber + * @brief Contains the @ref New_cell_constructor and @ref Pool_cell_constructor structures. + */ + #ifndef PM_COLUMN_CELL_CONSTRUCTORS_H #define PM_COLUMN_CELL_CONSTRUCTORS_H -#include //std::swap +#include //std::swap #include namespace Gudhi { namespace persistence_matrix { -template -struct New_cell_constructor{ - New_cell_constructor() {} - - template - Cell* construct(U&&...u) const { return new Cell(std::forward(u)...); } - - void destroy(Cell* cell) const { delete cell; } - - friend void swap(New_cell_constructor& col1, New_cell_constructor& col2){} +/** + * @brief @ref Cell factory. Constructs and destroyes 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 Destroyes 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) {} }; -template -struct Pool_cell_constructor{ -public: - Pool_cell_constructor() : cellPool_() {} - Pool_cell_constructor(const Pool_cell_constructor& col) : cellPool_(col.cellPool_) {} - Pool_cell_constructor(Pool_cell_constructor&& col) : cellPool_(std::move(col.cellPool_)) {} - - template - Cell* construct(U&&...u) { - return cellPool_.construct(std::forward(u)...); - } - - void destroy(Cell* cell) { - cellPool_.destroy(cell); - } - - Pool_cell_constructor& operator=(const Pool_cell_constructor& other){ - cellPool_ = other.cellPool_; - return *this; - } - - friend void swap(Pool_cell_constructor& col1, Pool_cell_constructor& col2){ - std::swap(col1.cellPool_, col2.cellPool_); - } - -private: - Simple_object_pool cellPool_; +/** + * @brief @ref Cell factory. Uses @ref Simple_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 Destroyes 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 +} // namespace persistence_matrix +} // namespace Gudhi -#endif // PM_COLUMN_CELL_CONSTRUCTORS_H +#endif // PM_COLUMN_CELL_CONSTRUCTORS_H diff --git a/src/Persistence_matrix/include/gudhi/Persistence_matrix/columns/cell_types.h b/src/Persistence_matrix/include/gudhi/Persistence_matrix/columns/cell_types.h index 489ff588b4..f184ad02ef 100644 --- a/src/Persistence_matrix/include/gudhi/Persistence_matrix/columns/cell_types.h +++ b/src/Persistence_matrix/include/gudhi/Persistence_matrix/columns/cell_types.h @@ -2,201 +2,322 @@ * 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-24 Inria * * Modification(s): * - YYYY/MM Author: Description of the modification */ +/** + * @file cell_types.h + * @author Hannah Schreiber + * @brief Contains the @ref Cell, @ref Cell_column_index and @ref Cell_field_element classes, as well as the + * @ref Dummy_cell_column_index_mixin and @ref Dummy_cell_field_element_mixin structures. Also defines the + * std::hash method for @ref Cell. + */ + #ifndef PM_MATRIX_CELL_H #define PM_MATRIX_CELL_H -#include //std::swap, std::exchange & std::move -#include //std::hash +#include //std::swap, std::exchange & std::move +#include //std::hash namespace Gudhi { namespace persistence_matrix { -struct Dummy_cell_column_index_mixin{ - Dummy_cell_column_index_mixin(){} - template - Dummy_cell_column_index_mixin([[maybe_unused]] index columnIndex){} - Dummy_cell_column_index_mixin([[maybe_unused]] const Dummy_cell_column_index_mixin& cell){}; - Dummy_cell_column_index_mixin([[maybe_unused]] Dummy_cell_column_index_mixin&& cell){}; +/** + * @brief Empty structure. + * Inheritated 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) {} + // Dummy_cell_column_index_mixin([[maybe_unused]] const Dummy_cell_column_index_mixin& cell){}; + // Dummy_cell_column_index_mixin([[maybe_unused]] Dummy_cell_column_index_mixin&& cell){}; - Dummy_cell_column_index_mixin& operator=([[maybe_unused]] const Dummy_cell_column_index_mixin& other){ - return *this; - }; + // Dummy_cell_column_index_mixin& operator=([[maybe_unused]] const Dummy_cell_column_index_mixin& other) { + // return *this; + // }; }; -struct Dummy_cell_field_element_mixin{ - Dummy_cell_field_element_mixin(){} - template - Dummy_cell_field_element_mixin([[maybe_unused]] Field_element_type t){} - Dummy_cell_field_element_mixin([[maybe_unused]] const Dummy_cell_field_element_mixin& cell){}; - Dummy_cell_field_element_mixin([[maybe_unused]] Dummy_cell_field_element_mixin&& cell){}; +/** + * @brief Empty structure. + * Inheritated 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_type t) {} + // Dummy_cell_field_element_mixin([[maybe_unused]] const Dummy_cell_field_element_mixin& cell){}; + // Dummy_cell_field_element_mixin([[maybe_unused]] Dummy_cell_field_element_mixin&& cell){}; - Dummy_cell_field_element_mixin& operator=([[maybe_unused]] const Dummy_cell_field_element_mixin& other){ - return *this; - }; + // Dummy_cell_field_element_mixin& operator=([[maybe_unused]] const Dummy_cell_field_element_mixin& other) { + // return *this; + // }; }; -template -class Cell_column_index +/** + * @brief Class managing the column index access of a cell. + * + * @tparam index MatIdx index type. + */ +template +class Cell_column_index { -public: - Cell_column_index() : columnIndex_(-1) {}; - Cell_column_index(index columnIndex) : columnIndex_(columnIndex) - {}; - Cell_column_index(const Cell_column_index& cell) : columnIndex_(cell.columnIndex_) - {}; - Cell_column_index(Cell_column_index&& cell) noexcept - : columnIndex_(std::exchange(cell.columnIndex_, 0)) - {}; - - index get_column_index() const{ - return columnIndex_; - }; - void set_column_index(index columnIndex){ - columnIndex_ = columnIndex; - } - - Cell_column_index& operator=(Cell_column_index other){ - std::swap(columnIndex_, other.columnIndex_); - return *this; - }; - -private: - index columnIndex_; + 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 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. */ }; -template -class Cell_field_element +/** + * @brief Class managing the value access of a cell. + * + * @tparam Field_element_type Type of a cell value. + */ +template +class Cell_field_element { -public: - Cell_field_element() : element_(0) {}; - Cell_field_element(Field_element_type element) : element_(element) - {}; - Cell_field_element(const Cell_field_element& cell) : element_(cell.element_) - {}; - Cell_field_element(Cell_field_element&& cell) noexcept - : element_(std::move(cell.element_)) - {}; - - Field_element_type& get_element(){ - return element_; - }; - const Field_element_type& get_element() const{ - return element_; - }; - void set_element(const Field_element_type &element){ - element_ = element; - } - - Cell_field_element& operator=(Cell_field_element other){ - std::swap(element_, other.element_); - return *this; - }; - -private: - Field_element_type element_; + public: + /** + * @brief Default constructor. Sets to the element to 0. + */ + Cell_field_element() : element_(0){}; + /** + * @brief Stores the given element. + * + * @param columnIndex Value to store. + */ + Cell_field_element(Field_element_type 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_type& get_element() { return element_; }; + /** + * @brief Returns the value stored in the cell. + * + * @return Const reference to the value of the cell. + */ + const Field_element_type& get_element() const { return element_; }; + /** + * @brief Sets the value. + * + * @param element Value to store. + */ + void set_element(const Field_element_type& element) { element_ = element; } + + /** + * @brief Assign operator. + */ + Cell_field_element& operator=(Cell_field_element other) { + std::swap(element_, other.element_); + return *this; + }; + + private: + Field_element_type element_; /**< Value of the cell. */ }; -template -class Cell : public Master_matrix::Cell_column_index_option, - public Master_matrix::Cell_field_element_option, - public Master_matrix::row_hook_type, - public Master_matrix::column_hook_type +/** + * @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 explicited in the matrix. + * + * @tparam Master_matrix An instanciation 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_type, + public Master_matrix::column_hook_type { -private: - using col_opt = typename Master_matrix::Cell_column_index_option; - using field_opt = typename Master_matrix::Cell_field_element_option; - -public: - using index = typename Master_matrix::index; - using id_index = typename Master_matrix::id_index; - using Field_element_type = typename Master_matrix::element_type; - - // using base_hook_matrix_row = typename Master_matrix::row_hook_type; //temporary during factorization process - - Cell(){}; - Cell(id_index rowIndex) - : col_opt(), - field_opt(), - rowIndex_(rowIndex) - {}; - Cell(index columnIndex, id_index rowIndex) - : col_opt(columnIndex), - field_opt(), - rowIndex_(rowIndex) - {}; - Cell(const Cell& cell) - : col_opt(static_cast(cell)), - field_opt(static_cast(cell)), - rowIndex_(cell.rowIndex_) - {}; - Cell(Cell&& cell) noexcept - : col_opt(std::move(static_cast(cell))), - field_opt(std::move(static_cast(cell))), - rowIndex_(std::exchange(cell.rowIndex_, 0)) - {}; - - id_index get_row_index() const{ - return rowIndex_; - }; - void set_row_index(id_index rowIndex){ - rowIndex_ = rowIndex; - }; - - Cell& operator=(Cell other){ - col_opt::operator=(other); - field_opt::operator=(other); - std::swap(rowIndex_, other.rowIndex_); - return *this; - }; - - friend bool operator<(const Cell& c1, const Cell& c2){ - return c1.get_row_index() < c2.get_row_index(); - } - friend bool operator==(const Cell& c1, const Cell& c2){ - return c1.get_row_index() == c2.get_row_index(); - } - - // operator unsigned int() const{ - // return rowIndex_; - // } - // operator std::pair() const{ - // if constexpr (Master_matrix::Option_list::is_z2){ - // return {rowIndex_, 1}; - // } else { - // return {rowIndex_, field_opt::element_}; - // } - // } - operator id_index() const{ - return rowIndex_; - } - 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_; + private: + using col_opt = typename Master_matrix::Cell_column_index_option; + using field_opt = typename Master_matrix::Cell_field_element_option; + + public: + using index = typename Master_matrix::index; /**< Column index type. */ + using id_index = typename Master_matrix::id_index; /**< Row index type. */ + using Field_element_type = typename Master_matrix::element_type; /**< 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 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 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 Row index of the cell. + */ + id_index get_row_index() const { return rowIndex_; }; + /** + * @brief Sets the row index stored in the cell. + * + * @param 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 +} // namespace persistence_matrix +} // namespace Gudhi -template -struct std::hash > -{ - size_t operator()(const Gudhi::persistence_matrix::Cell& cell) const - { - return std::hash()(cell.get_row_index()); - } +/** + * @brief Hash method for @ref Gudhi::persistence_matrix::Cell. + * + * The cells are differentiated by their row indices only. For exemple, 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 > { + size_t operator()(const Gudhi::persistence_matrix::Cell& cell) const { + return std::hash()(cell.get_row_index()); + } }; -#endif // PM_MATRIX_CELL_H +#endif // PM_MATRIX_CELL_H diff --git a/src/Persistence_matrix/include/gudhi/Persistence_matrix/columns/chain_column_extra_properties.h b/src/Persistence_matrix/include/gudhi/Persistence_matrix/columns/chain_column_extra_properties.h index d8e3de9809..f0901dfc67 100644 --- a/src/Persistence_matrix/include/gudhi/Persistence_matrix/columns/chain_column_extra_properties.h +++ b/src/Persistence_matrix/include/gudhi/Persistence_matrix/columns/chain_column_extra_properties.h @@ -2,72 +2,145 @@ * 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-24 Inria * * Modification(s): * - YYYY/MM Author: Description of the modification */ +/** + * @file chain_column_extra_properties.h + * @author Hannah Schreiber + * @brief Contains the @ref Chain_column_extra_properties class and @ref Dummy_chain_properties structure. + */ + #ifndef PM_CHAIN_COLUMN_PROP_H #define PM_CHAIN_COLUMN_PROP_H -#include //std::swap +#include //std::swap namespace Gudhi { namespace persistence_matrix { -struct Dummy_chain_properties{ - Dummy_chain_properties(){} - Dummy_chain_properties([[maybe_unused]] int pivot){} - Dummy_chain_properties([[maybe_unused]] int pivot, [[maybe_unused]] int pair){} - Dummy_chain_properties([[maybe_unused]] const Dummy_chain_properties& col){} - Dummy_chain_properties([[maybe_unused]] Dummy_chain_properties&& col){} +/** + * @brief Empty structure. + * Inheritated instead of @ref Chain_column_extra_properties, when the columns are not meant for chain matrices. + */ +struct Dummy_chain_properties +{ + // Dummy_chain_properties() {} + // Dummy_chain_properties([[maybe_unused]] int pivot) {} + Dummy_chain_properties([[maybe_unused]] int pivot = 0, [[maybe_unused]] int pair = 0) {} + // Dummy_chain_properties([[maybe_unused]] const Dummy_chain_properties& col) {} + // Dummy_chain_properties([[maybe_unused]] Dummy_chain_properties&& col) {} - Dummy_chain_properties& operator=([[maybe_unused]] const Dummy_chain_properties& other){ return *this; } + // Dummy_chain_properties& operator=([[maybe_unused]] const Dummy_chain_properties& other) { return *this; } - friend void swap([[maybe_unused]] Dummy_chain_properties& col1, [[maybe_unused]] Dummy_chain_properties& col2){} + friend void swap([[maybe_unused]] Dummy_chain_properties& col1, [[maybe_unused]] Dummy_chain_properties& col2) {} }; -template -class Chain_column_extra_properties{ -public: - using index = typename Master_matrix::index; - using id_index = typename Master_matrix::id_index; - - Chain_column_extra_properties() : pivot_(-1), pairedColumn_(-1) {} - Chain_column_extra_properties(id_index pivot) : pivot_(pivot), pairedColumn_(-1) {} - Chain_column_extra_properties(id_index pivot, index pair) : pivot_(pivot), pairedColumn_(pair) {} - Chain_column_extra_properties(const Chain_column_extra_properties& col) - : pivot_(col.pivot_), pairedColumn_(col.pairedColumn_) {} - Chain_column_extra_properties(Chain_column_extra_properties&& col) - : pivot_(std::exchange(col.pivot_, -1)), pairedColumn_(std::exchange(col.pairedColumn_, -1)) {} +/** + * @brief Class managing the pivot and partitioning of columns in @ref Chain_matrix. + * + * The columns of a chain matrix are partitioned in three sets: \f$ F \f$, \f$ G \f$ and \f$ H \f$ with a + * bijection between \f$ G \f$ and \f$ H \f$. If a column is in \f$ F \f$, the value of + * @ref Chain_column_extra_properties::pairedColumn_ is -1, while the value corresponds to the MatIdx index of + * the image of the bijection if the column is in either \f$ G \f$ or \f$ H \f$. See [TODO: zigzag paper] for + * more details. + * + * @tparam Master_matrix An instanciation of @ref Matrix from which all types and options are deduced. + */ +template +class Chain_column_extra_properties +{ + public: + using index = typename Master_matrix::index; /**< MatIdx index type. */ + using id_index = typename Master_matrix::id_index; /**< IDIdx index type. */ - index get_paired_chain_index() const { return pairedColumn_; } - bool is_paired() const { return pairedColumn_ != static_cast(-1); } - void assign_paired_chain(index other_col){ pairedColumn_ = other_col; } - void unassign_paired_chain() { pairedColumn_ = -1; }; + /** + * @brief Default constructor. Sets the pivot and pair to -1, which means "not existing". + */ + Chain_column_extra_properties() : pivot_(-1), pairedColumn_(-1) {} + /** + * @brief Constructor setting the pivot at the given value and the pair to -1 (i.e. not paired). + * + * @param pivot Row index of the pivot. Corresponds to the IDIdx index of the face represented by the column. + */ + Chain_column_extra_properties(id_index pivot) : pivot_(pivot), pairedColumn_(-1) {} + /** + * @brief Constructor setting the pivot and the pair at the given values. + * + * @param pivot Row index of the pivot. Corresponds to the IDIdx index of the face represented by the column. + * @param pair MatIdx index of the pair of the column. + */ + Chain_column_extra_properties(id_index pivot, index pair) : pivot_(pivot), pairedColumn_(pair) {} + /** + * @brief Copy constructor. + * + * @param col Column to copy. + */ + Chain_column_extra_properties(const Chain_column_extra_properties& col) + : pivot_(col.pivot_), pairedColumn_(col.pairedColumn_) {} + /** + * @brief Move constructor. + * + * @param col Column to move. + */ + Chain_column_extra_properties(Chain_column_extra_properties&& col) + : pivot_(std::exchange(col.pivot_, -1)), pairedColumn_(std::exchange(col.pairedColumn_, -1)) {} - Chain_column_extra_properties& operator=(const Chain_column_extra_properties& other){ - pivot_ = other.pivot_; - pairedColumn_ = other.pairedColumn_; - return *this; - } + /** + * @brief Returns -1 if the column is not paired, the matIdx of the pair otherwise. + * + * @return -1 if the column is not paired, the matIdx of the pair otherwise. + */ + index get_paired_chain_index() const { return pairedColumn_; } + /** + * @brief Indicates if the column is paired or not. + * + * @return true If the column is paired. + * @return false Otherwise. + */ + bool is_paired() const { return pairedColumn_ != static_cast(-1); } + /** + * @brief Sets the value of the pair. + * + * @param other_col MatIdx of the pair column. + */ + void assign_paired_chain(index other_col) { pairedColumn_ = other_col; } + /** + * @brief Unpairs a column. + */ + void unassign_paired_chain() { pairedColumn_ = -1; }; - friend void swap(Chain_column_extra_properties& col1, Chain_column_extra_properties& col2){ - std::swap(col1.pivot_, col2.pivot_); - std::swap(col1.pairedColumn_, col2.pairedColumn_); - } + /** + * @brief Assign operator. + */ + Chain_column_extra_properties& operator=(const Chain_column_extra_properties& other) { + pivot_ = other.pivot_; + pairedColumn_ = other.pairedColumn_; + return *this; + } + /** + * @brief Swap operator. + */ + friend void swap(Chain_column_extra_properties& col1, Chain_column_extra_properties& col2) { + std::swap(col1.pivot_, col2.pivot_); + std::swap(col1.pairedColumn_, col2.pairedColumn_); + } -protected: - id_index get_pivot() const { return pivot_; } - void swap_pivots(Chain_column_extra_properties& other) { std::swap(pivot_, other.pivot_); } + protected: + id_index get_pivot() const { return pivot_; } + void swap_pivots(Chain_column_extra_properties& other) { std::swap(pivot_, other.pivot_); } -private: - id_index pivot_; //simplex index associated to the chain - index pairedColumn_; //represents the (F, G x H) partition of the columns + private: + id_index pivot_; /**< IDIdx index associated to the chain */ + index pairedColumn_; /**< Represents the (F, G x H) partition of the columns. + -1 if in F, MatIdx of image of bijection otherwise. + The pivot of a column in G will always be smaller than the pivot of its image in H. */ }; -} //namespace persistence_matrix -} //namespace Gudhi +} // namespace persistence_matrix +} // namespace Gudhi -#endif // PM_CHAIN_COLUMN_PROP_H +#endif // PM_CHAIN_COLUMN_PROP_H diff --git a/src/Persistence_matrix/include/gudhi/Persistence_matrix/columns/column_dimension_holder.h b/src/Persistence_matrix/include/gudhi/Persistence_matrix/columns/column_dimension_holder.h index 45fb6ed18d..2aa7febef7 100644 --- a/src/Persistence_matrix/include/gudhi/Persistence_matrix/columns/column_dimension_holder.h +++ b/src/Persistence_matrix/include/gudhi/Persistence_matrix/columns/column_dimension_holder.h @@ -2,60 +2,103 @@ * 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-24 Inria * * Modification(s): * - YYYY/MM Author: Description of the modification */ +/** + * @file column_dimension_holder.h + * @author Hannah Schreiber + * @brief Contains the @ref Column_dimension_holder class and @ref Dummy_dimension_holder structure. + */ + #ifndef PM_COLUMN_DIM_HOLDER_H #define PM_COLUMN_DIM_HOLDER_H -#include //std::swap +#include //std::swap namespace Gudhi { namespace persistence_matrix { -struct Dummy_dimension_holder{ - Dummy_dimension_holder(){} - template - Dummy_dimension_holder([[maybe_unused]] dimension_type dim){} - Dummy_dimension_holder([[maybe_unused]] const Dummy_dimension_holder& col){} - Dummy_dimension_holder([[maybe_unused]] Dummy_dimension_holder&& col){} +/** + * @brief Empty structure. + * Inheritated instead of @ref Column_dimension_holder, when the columns are not storing a dimension. + */ +struct Dummy_dimension_holder +{ + Dummy_dimension_holder() {} + template + Dummy_dimension_holder([[maybe_unused]] dimension_type dim) {} + // Dummy_dimension_holder([[maybe_unused]] const Dummy_dimension_holder& col) {} + // Dummy_dimension_holder([[maybe_unused]] Dummy_dimension_holder&& col) {} - Dummy_dimension_holder& operator=([[maybe_unused]] const Dummy_dimension_holder& other){ return *this; } + // Dummy_dimension_holder& operator=([[maybe_unused]] const Dummy_dimension_holder& other) { return *this; } - friend void swap([[maybe_unused]] Dummy_dimension_holder& col1, [[maybe_unused]] Dummy_dimension_holder& col2){} + friend void swap([[maybe_unused]] Dummy_dimension_holder& col1, [[maybe_unused]] Dummy_dimension_holder& col2) {} }; -template -struct Column_dimension_holder{ - using dimension_type = typename Master_matrix::dimension_type; - - Column_dimension_holder() : dim_(Master_matrix::Option_list::is_of_boundary_type ? 0 : -1) {} - Column_dimension_holder(dimension_type dim) : dim_(dim) {} - Column_dimension_holder(const Column_dimension_holder& col) : dim_(col.dim_) {} - Column_dimension_holder(Column_dimension_holder&& col) : dim_(std::exchange(col.dim_, -1)) {} +/** + * @brief Class managing the dimension access of a column. + * + * @tparam Master_matrix An instanciation of @ref Matrix from which all types and options are deduced. + */ +template +struct Column_dimension_holder +{ + using dimension_type = typename Master_matrix::dimension_type; /**< Dimension value type. */ - dimension_type get_dimension() const { return dim_; } + /** + * @brief Default constructor. Sets the dimension to 0 for boundary matrices and to -1 for chain matrices. + */ + Column_dimension_holder() : dim_(Master_matrix::Option_list::is_of_boundary_type ? 0 : -1) {} + /** + * @brief Constructor setting the dimension to the given value. + * + * @param dim Dimension of the column. + */ + Column_dimension_holder(dimension_type dim) : dim_(dim) {} + /** + * @brief Copy constructor. + * + * @param col Column to copy. + */ + Column_dimension_holder(const Column_dimension_holder& col) : dim_(col.dim_) {} + /** + * @brief Move constructor. + * + * @param col Column to move. + */ + Column_dimension_holder(Column_dimension_holder&& col) : dim_(std::exchange(col.dim_, -1)) {} - Column_dimension_holder& operator=(const Column_dimension_holder& other){ - dim_ = other.dim_; - return *this; - } + /** + * @brief Returns the dimension of the column. + * + * @return The dimension of the column. + */ + dimension_type get_dimension() const { return dim_; } - friend void swap(Column_dimension_holder& col1, Column_dimension_holder& col2){ - std::swap(col1.dim_, col2.dim_); - } + /** + * @brief Assign operator. + */ + Column_dimension_holder& operator=(const Column_dimension_holder& other) { + dim_ = other.dim_; + return *this; + } + /** + * @brief Swap operator. + */ + friend void swap(Column_dimension_holder& col1, Column_dimension_holder& col2) { std::swap(col1.dim_, col2.dim_); } -protected: - void swap_dimension(Column_dimension_holder& other) { std::swap(dim_, other.dim_); } + protected: + void swap_dimension(Column_dimension_holder& other) { std::swap(dim_, other.dim_); } -private: - dimension_type dim_; + private: + dimension_type dim_; /**< Dimension of the column. */ }; -} //namespace persistence_matrix -} //namespace Gudhi +} // namespace persistence_matrix +} // namespace Gudhi -#endif // PM_COLUMN_DIM_HOLDER_H +#endif // PM_COLUMN_DIM_HOLDER_H diff --git a/src/Persistence_matrix/include/gudhi/Persistence_matrix/columns/heap_column.h b/src/Persistence_matrix/include/gudhi/Persistence_matrix/columns/heap_column.h index efa06e396c..84a2677917 100644 --- a/src/Persistence_matrix/include/gudhi/Persistence_matrix/columns/heap_column.h +++ b/src/Persistence_matrix/include/gudhi/Persistence_matrix/columns/heap_column.h @@ -2,20 +2,26 @@ * 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-24 Inria * * Modification(s): * - YYYY/MM Author: Description of the modification */ +/** + * @file heap_column.h + * @author Hannah Schreiber + * @brief Contains the @ref Heap_column class. Also defines the std::hash method for @ref Heap_column. + */ + #ifndef PM_HEAP_COLUMN_H #define PM_HEAP_COLUMN_H #include #include #include -#include //binary_search -#include //std::swap, std::move & std::exchange +#include //binary_search +#include //std::swap, std::move & std::exchange #include @@ -24,984 +30,1041 @@ namespace Gudhi { namespace persistence_matrix { -template > -class Heap_column : public Master_matrix::Column_dimension_option, - public Master_matrix::Chain_column_option -{ -public: - using Master = Master_matrix; - using Field_operators = typename Master_matrix::Field_operators; - using Field_element_type = typename Master_matrix::element_type; - using index = typename Master_matrix::index; - using id_index = typename Master_matrix::id_index; - using dimension_type = typename Master_matrix::dimension_type; - - using Cell = typename Master_matrix::Cell_type; - using Column_type = std::vector; - using iterator = boost::indirect_iterator; - using const_iterator = boost::indirect_iterator; - using reverse_iterator = boost::indirect_iterator; - using const_reverse_iterator = boost::indirect_iterator; - - Heap_column(Field_operators* operators = nullptr, Cell_constructor* cellConstructor = nullptr); - template - Heap_column(const Container_type& nonZeroRowIndices, Field_operators* operators, Cell_constructor* cellConstructor); //has to be a boundary for boundary, has no sense for chain if dimension is needed - template - Heap_column(const Container_type& nonZeroChainRowIndices, dimension_type dimension, Field_operators* operators, Cell_constructor* cellConstructor); //dimension gets ignored for base - Heap_column(const Heap_column& column, Field_operators* operators = nullptr, Cell_constructor* cellConstructor = nullptr); - Heap_column(Heap_column&& column) noexcept; - ~Heap_column(); - - //just for the sake of the interface - //row containers and column index are ignored as row access is not implemented for heap columns - template - Heap_column(index columnIndex, const Container_type& nonZeroRowIndices, Row_container_type *rowContainer, Field_operators* operators, Cell_constructor* cellConstructor); - template - Heap_column(index columnIndex, const Container_type& nonZeroChainRowIndices, dimension_type dimension, Row_container_type *rowContainer, Field_operators* operators, Cell_constructor* cellConstructor); - template - Heap_column(const Heap_column& column, index columnIndex, Row_container_type *rowContainer, Field_operators* operators = nullptr, Cell_constructor* cellConstructor = nullptr); - - std::vector get_content(int columnLength = -1) const; - bool is_non_zero(id_index rowIndex) const; - bool is_empty(); - std::size_t size() const; //size of the underlying vector, not the number of non zero row cells. - - //**************** - //only for base and boundary - template - void reorder(const Map_type& valueMap, [[maybe_unused]] index columnIndex = -1); //used for lazy row swaps - void clear(); - void clear(id_index rowIndex); - //**************** - - //**************** - //only for chain and boundary - id_index get_pivot(); - Field_element_type get_pivot_value(); - //**************** - - iterator begin() noexcept; - const_iterator begin() const noexcept; - iterator end() noexcept; - const_iterator end() const noexcept; - reverse_iterator rbegin() noexcept; - const_reverse_iterator rbegin() const noexcept; - reverse_iterator rend() noexcept; - const_reverse_iterator rend() const noexcept; - - template - Heap_column& operator+=(const Cell_range& column); //for base & boundary except vector //may not work if Cell_range = Heap_column - Heap_column& operator+=(Heap_column &column); //for chain and vector - - Heap_column& operator*=(unsigned int v); - - //this = v * this + column - template - Heap_column& multiply_and_add(const Field_element_type& val, const Cell_range& column); //for base & boundary except vector //may not work if Cell_range = Heap_column - Heap_column& multiply_and_add(const Field_element_type& val, Heap_column& column); //for chain and vector - //this = this + column * v - template - Heap_column& multiply_and_add(const Cell_range& column, const Field_element_type& val); //for base & boundary except vector //may not work if Cell_range = Heap_column - Heap_column& multiply_and_add(Heap_column& column, const Field_element_type& val); //for chain and vector - - friend bool operator==(const Heap_column& c1, const Heap_column& c2){ - if (&c1 == &c2) return true; - return c1.get_content() == c2.get_content(); - } - friend bool operator<(const Heap_column& c1, const Heap_column& c2){ - if (&c1 == &c2) return false; - - auto cont1 = c1.get_content(); - auto cont2 = c2.get_content(); - - auto it1 = cont1.begin(); - auto it2 = cont2.begin(); - while (it1 != cont1.end() && it2 != cont2.end()) { - if (*it1 == 0u && *it2 != 0u) return false; - if (*it1 != 0u && *it2 == 0u) return true; - if (*it1 != *it2) return *it1 < *it2; - ++it1; ++it2; - } - return it2 != cont2.end(); - } - - // void set_operators(Field_operators* operators){ operators_ = operators; } - - //Disabled with row access. - Heap_column& operator=(const Heap_column& other); - - friend void swap(Heap_column& col1, Heap_column& col2){ - swap(static_cast(col1), - static_cast(col2)); - swap(static_cast(col1), - static_cast(col2)); - col1.column_.swap(col2.column_); - std::swap(col1.insertsSinceLastPrune_, col2.insertsSinceLastPrune_); - std::swap(col1.operators_, col2.operators_); - std::swap(col1.cellPool_, col2.cellPool_); - } - -private: - using dim_opt = typename Master_matrix::Column_dimension_option; - using chain_opt = typename Master_matrix::Chain_column_option; - - struct{ - bool operator()(const Cell* c1, const Cell* c2) const { return *c1 < *c2; } - } cellPointerComp_; - - Column_type column_; - unsigned int insertsSinceLastPrune_; - Field_operators* operators_; - Cell_constructor* cellPool_; - - void _prune(); - Cell* _pop_pivot(); - template - bool _add(const Cell_range& column); - template - bool _multiply_and_add(const Field_element_type& val, const Cell_range& column); - template - bool _multiply_and_add(const Cell_range& column, const Field_element_type& val); - - void _verifyCellConstructor(){ - if (cellPool_ == nullptr){ - if constexpr (std::is_same_v >){ - cellPool_ = &Master_matrix::defaultCellConstructor; - } else { - throw std::invalid_argument("Cell constructor pointer cannot be null."); - } - } - } +/** + * @brief Column class following the @ref PersistenceMatrixColumn concept. Not compatible with row access. + * + * Column based on a heap structure. The heap is represented as a vector sorted as a heap. The top of the heap is + * the cell with the biggest row index. The sum of two columns is lazy: the content of the source is simply inserted + * into the heap of the target. Therefore the underlying vector can contain several cells with the same row index. + * The real value of a cell at a row index corresponds to the sum in the coeffcient field of all values with same + * row index. + * + * @tparam Master_matrix An instanciation of @ref Matrix from which all types and options are deduced. + * @tparam Cell_constructor Factory of @ref Cell classes. + */ +template > +class Heap_column : public Master_matrix::Column_dimension_option, public Master_matrix::Chain_column_option +{ + public: + using Master = Master_matrix; + using Field_operators = typename Master_matrix::Field_operators; + using Field_element_type = typename Master_matrix::element_type; + using index = typename Master_matrix::index; + using id_index = typename Master_matrix::id_index; + using dimension_type = typename Master_matrix::dimension_type; + using Cell = typename Master_matrix::Cell_type; + using Column_type = std::vector; + using iterator = boost::indirect_iterator; + using const_iterator = boost::indirect_iterator; + using reverse_iterator = boost::indirect_iterator; + using const_reverse_iterator = boost::indirect_iterator; + + Heap_column(Field_operators* operators = nullptr, Cell_constructor* cellConstructor = nullptr); + template + Heap_column(const Container_type& nonZeroRowIndices, Field_operators* operators, Cell_constructor* cellConstructor); + template + Heap_column(const Container_type& nonZeroChainRowIndices, + dimension_type dimension, + Field_operators* operators, + Cell_constructor* cellConstructor); + Heap_column(const Heap_column& column, Field_operators* operators = nullptr, + Cell_constructor* cellConstructor = nullptr); + Heap_column(Heap_column&& column) noexcept; + ~Heap_column(); + + // just for the sake of the interface + // row containers and column index are ignored as row access is not implemented for heap columns + template + Heap_column(index columnIndex, + const Container_type& nonZeroRowIndices, + Row_container_type* rowContainer, + Field_operators* operators, + Cell_constructor* cellConstructor); + template + Heap_column(index columnIndex, + const Container_type& nonZeroChainRowIndices, + dimension_type dimension, + Row_container_type* rowContainer, + Field_operators* operators, + Cell_constructor* cellConstructor); + template + Heap_column(const Heap_column& column, + index columnIndex, + Row_container_type* rowContainer, + Field_operators* operators = nullptr, + Cell_constructor* cellConstructor = nullptr); + + std::vector get_content(int columnLength = -1) const; + bool is_non_zero(id_index rowIndex) const; + bool is_empty(); + std::size_t size() const; + + template + void reorder(const Map_type& valueMap, [[maybe_unused]] index columnIndex = -1); + void clear(); + void clear(id_index rowIndex); + + id_index get_pivot(); + Field_element_type get_pivot_value(); + + iterator begin() noexcept; + const_iterator begin() const noexcept; + iterator end() noexcept; + const_iterator end() const noexcept; + reverse_iterator rbegin() noexcept; + const_reverse_iterator rbegin() const noexcept; + reverse_iterator rend() noexcept; + const_reverse_iterator rend() const noexcept; + + template + Heap_column& operator+=(const Cell_range& column); + Heap_column& operator+=(Heap_column& column); + + Heap_column& operator*=(unsigned int v); + + // this = v * this + column + template + Heap_column& multiply_and_add(const Field_element_type& val, const Cell_range& column); + Heap_column& multiply_and_add(const Field_element_type& val, Heap_column& column); + // this = this + column * v + template + Heap_column& multiply_and_add(const Cell_range& column, const Field_element_type& val); + Heap_column& multiply_and_add(Heap_column& column, const Field_element_type& val); + + friend bool operator==(const Heap_column& c1, const Heap_column& c2) { + if (&c1 == &c2) return true; + return c1.get_content() == c2.get_content(); + } + friend bool operator<(const Heap_column& c1, const Heap_column& c2) { + if (&c1 == &c2) return false; + + auto cont1 = c1.get_content(); + auto cont2 = c2.get_content(); + + auto it1 = cont1.begin(); + auto it2 = cont2.begin(); + while (it1 != cont1.end() && it2 != cont2.end()) { + if (*it1 == 0u && *it2 != 0u) return false; + if (*it1 != 0u && *it2 == 0u) return true; + if (*it1 != *it2) return *it1 < *it2; + ++it1; + ++it2; + } + return it2 != cont2.end(); + } + + // void set_operators(Field_operators* operators){ operators_ = operators; } + + // Disabled with row access. + Heap_column& operator=(const Heap_column& other); + + friend void swap(Heap_column& col1, Heap_column& col2) { + swap(static_cast(col1), + static_cast(col2)); + swap(static_cast(col1), + static_cast(col2)); + col1.column_.swap(col2.column_); + std::swap(col1.insertsSinceLastPrune_, col2.insertsSinceLastPrune_); + std::swap(col1.operators_, col2.operators_); + std::swap(col1.cellPool_, col2.cellPool_); + } + + private: + using dim_opt = typename Master_matrix::Column_dimension_option; + using chain_opt = typename Master_matrix::Chain_column_option; + + struct { + bool operator()(const Cell* c1, const Cell* c2) const { return *c1 < *c2; } + } cellPointerComp_; + + Column_type column_; + unsigned int insertsSinceLastPrune_; + Field_operators* operators_; + Cell_constructor* cellPool_; + + void _prune(); + Cell* _pop_pivot(); + template + bool _add(const Cell_range& column); + template + bool _multiply_and_add(const Field_element_type& val, const Cell_range& column); + template + bool _multiply_and_add(const Cell_range& column, const Field_element_type& val); + + void _verifyCellConstructor() { + if (cellPool_ == nullptr) { + if constexpr (std::is_same_v >) { + cellPool_ = &Master_matrix::defaultCellConstructor; + } else { + throw std::invalid_argument("Cell constructor pointer cannot be null."); + } + } + } }; -template -inline Heap_column::Heap_column(Field_operators* operators, Cell_constructor* cellConstructor) - : dim_opt(), chain_opt(), insertsSinceLastPrune_(0), operators_(operators) -{ - if (operators_ == nullptr && cellPool_ == nullptr) return; //to allow default constructor which gives a dummy column - _verifyCellConstructor(); -} - -template -template -inline Heap_column::Heap_column(const Container_type &nonZeroRowIndices, Field_operators* operators, Cell_constructor* cellConstructor) - : dim_opt(nonZeroRowIndices.size() == 0 ? 0 : nonZeroRowIndices.size() - 1), - chain_opt(), - column_(nonZeroRowIndices.size(), nullptr), - insertsSinceLastPrune_(0), - operators_(operators), - cellPool_(cellConstructor) -{ - static_assert(!Master_matrix::isNonBasic || Master_matrix::Option_list::is_of_boundary_type, - "Constructor not available for chain columns, please specify the dimension of the chain."); - - _verifyCellConstructor(); - - index i = 0; - if constexpr (Master_matrix::Option_list::is_z2){ - for (id_index id : nonZeroRowIndices){ - column_[i++] = cellPool_->construct(id); - } - } else { - for (const auto& p : nonZeroRowIndices){ - column_[i] = cellPool_->construct(p.first); - column_[i++]->set_element(operators_->get_value(p.second)); - } - } - std::make_heap(column_.begin(), column_.end(), cellPointerComp_); -} - -template -template -inline Heap_column::Heap_column( - const Container_type &nonZeroRowIndices, dimension_type dimension, Field_operators* operators, Cell_constructor* cellConstructor) - : dim_opt(dimension), - chain_opt([&]{ - if constexpr (Master_matrix::Option_list::is_z2){ - return nonZeroRowIndices.begin() == nonZeroRowIndices.end() ? -1 : *std::prev(nonZeroRowIndices.end()); - } else { - return nonZeroRowIndices.begin() == nonZeroRowIndices.end() ? -1 : std::prev(nonZeroRowIndices.end())->first; - } - }()), - column_(nonZeroRowIndices.size(), nullptr), - insertsSinceLastPrune_(0), - operators_(operators), - cellPool_(cellConstructor) -{ - _verifyCellConstructor(); - - index i = 0; - if constexpr (Master_matrix::Option_list::is_z2){ - for (id_index id : nonZeroRowIndices){ - column_[i++] = cellPool_->construct(id); - } - } else { - for (const auto& p : nonZeroRowIndices){ - column_[i] = cellPool_->construct(p.first); - column_[i++]->set_element(operators_->get_value(p.second)); - } - } - std::make_heap(column_.begin(), column_.end(), cellPointerComp_); -} - -template -inline Heap_column::Heap_column(const Heap_column &column, Field_operators* operators, Cell_constructor* cellConstructor) - : dim_opt(static_cast(column)), - chain_opt(static_cast(column)), - column_(column.column_.size(), nullptr), - insertsSinceLastPrune_(0), - operators_(operators == nullptr ? column.operators_ : operators), - cellPool_(cellConstructor == nullptr ? column.cellPool_ : cellConstructor) -{ - static_assert(!Master_matrix::Option_list::has_row_access, - "Simple copy constructor not available when row access option enabled. Please specify the new column index and the row container."); - - index i = 0; - for (const Cell* cell : column.column_){ - if constexpr (Master_matrix::Option_list::is_z2){ - column_[i++] = cellPool_->construct(cell->get_row_index()); - } else { - column_[i] = cellPool_->construct(cell->get_row_index()); - column_[i++]->set_element(cell->get_element()); - } - } - //column.column_ already ordered as a heap, so no need of make_heap. -} - -template -inline Heap_column::Heap_column(Heap_column &&column) noexcept - : dim_opt(std::move(static_cast(column))), - chain_opt(std::move(static_cast(column))), - column_(std::move(column.column_)), - insertsSinceLastPrune_(std::exchange(column.insertsSinceLastPrune_, 0)), - operators_(std::exchange(column.operators_, nullptr)), - cellPool_(std::exchange(column.cellPool_, nullptr)) +template +inline Heap_column::Heap_column(Field_operators* operators, + Cell_constructor* cellConstructor) + : dim_opt(), chain_opt(), insertsSinceLastPrune_(0), operators_(operators) +{ + if (operators_ == nullptr && cellPool_ == nullptr) return; //to allow default constructor which gives a dummy column + _verifyCellConstructor(); +} + +template +template +inline Heap_column::Heap_column(const Container_type& nonZeroRowIndices, + Field_operators* operators, + Cell_constructor* cellConstructor) + : dim_opt(nonZeroRowIndices.size() == 0 ? 0 : nonZeroRowIndices.size() - 1), + chain_opt(), + column_(nonZeroRowIndices.size(), nullptr), + insertsSinceLastPrune_(0), + operators_(operators), + cellPool_(cellConstructor) +{ + static_assert(!Master_matrix::isNonBasic || Master_matrix::Option_list::is_of_boundary_type, + "Constructor not available for chain columns, please specify the dimension of the chain."); + + _verifyCellConstructor(); + + index i = 0; + if constexpr (Master_matrix::Option_list::is_z2) { + for (id_index id : nonZeroRowIndices) { + column_[i++] = cellPool_->construct(id); + } + } else { + for (const auto& p : nonZeroRowIndices) { + column_[i] = cellPool_->construct(p.first); + column_[i++]->set_element(operators_->get_value(p.second)); + } + } + std::make_heap(column_.begin(), column_.end(), cellPointerComp_); +} + +template +template +inline Heap_column::Heap_column(const Container_type& nonZeroRowIndices, + dimension_type dimension, + Field_operators* operators, + Cell_constructor* cellConstructor) + : dim_opt(dimension), + chain_opt([&] { + if constexpr (Master_matrix::Option_list::is_z2) { + return nonZeroRowIndices.begin() == nonZeroRowIndices.end() ? -1 : *std::prev(nonZeroRowIndices.end()); + } else { + return nonZeroRowIndices.begin() == nonZeroRowIndices.end() ? -1 : std::prev(nonZeroRowIndices.end())->first; + } + }()), + column_(nonZeroRowIndices.size(), nullptr), + insertsSinceLastPrune_(0), + operators_(operators), + cellPool_(cellConstructor) +{ + _verifyCellConstructor(); + + index i = 0; + if constexpr (Master_matrix::Option_list::is_z2) { + for (id_index id : nonZeroRowIndices) { + column_[i++] = cellPool_->construct(id); + } + } else { + for (const auto& p : nonZeroRowIndices) { + column_[i] = cellPool_->construct(p.first); + column_[i++]->set_element(operators_->get_value(p.second)); + } + } + std::make_heap(column_.begin(), column_.end(), cellPointerComp_); +} + +template +inline Heap_column::Heap_column(const Heap_column& column, + Field_operators* operators, + Cell_constructor* cellConstructor) + : dim_opt(static_cast(column)), + chain_opt(static_cast(column)), + column_(column.column_.size(), nullptr), + insertsSinceLastPrune_(0), + operators_(operators == nullptr ? column.operators_ : operators), + cellPool_(cellConstructor == nullptr ? column.cellPool_ : cellConstructor) +{ + static_assert(!Master_matrix::Option_list::has_row_access, + "Simple copy constructor not available when row access option enabled. Please specify the new column " + "index and the row container."); + + index i = 0; + for (const Cell* cell : column.column_) { + if constexpr (Master_matrix::Option_list::is_z2) { + column_[i++] = cellPool_->construct(cell->get_row_index()); + } else { + column_[i] = cellPool_->construct(cell->get_row_index()); + column_[i++]->set_element(cell->get_element()); + } + } + // column.column_ already ordered as a heap, so no need of make_heap. +} + +template +inline Heap_column::Heap_column(Heap_column&& column) noexcept + : dim_opt(std::move(static_cast(column))), + chain_opt(std::move(static_cast(column))), + column_(std::move(column.column_)), + insertsSinceLastPrune_(std::exchange(column.insertsSinceLastPrune_, 0)), + operators_(std::exchange(column.operators_, nullptr)), + cellPool_(std::exchange(column.cellPool_, nullptr)) {} -template -template -inline Heap_column::Heap_column( - index columnIndex, const Container_type &nonZeroRowIndices, Row_container_type *rowContainer, Field_operators* operators, Cell_constructor* cellConstructor) - : dim_opt(nonZeroRowIndices.size() == 0 ? 0 : nonZeroRowIndices.size() - 1), - chain_opt([&]{ - if constexpr (Master_matrix::Option_list::is_z2){ - return nonZeroRowIndices.begin() == nonZeroRowIndices.end() ? -1 : *std::prev(nonZeroRowIndices.end()); - } else { - return nonZeroRowIndices.begin() == nonZeroRowIndices.end() ? -1 : std::prev(nonZeroRowIndices.end())->first; - } - }()), - column_(nonZeroRowIndices.size(), nullptr), - insertsSinceLastPrune_(0), - operators_(operators), - cellPool_(cellConstructor) -{ - static_assert(!Master_matrix::isNonBasic || Master_matrix::Option_list::is_of_boundary_type, - "Constructor not available for chain columns, please specify the dimension of the chain."); - - _verifyCellConstructor(); - - index i = 0; - if constexpr (Master_matrix::Option_list::is_z2){ - for (id_index id : nonZeroRowIndices){ - column_[i++] = cellPool_->construct(id); - } - } else { - for (const auto& p : nonZeroRowIndices){ - column_[i] = cellPool_->construct(p.first); - column_[i++]->set_element(operators_->get_value(p.second)); - } - } - std::make_heap(column_.begin(), column_.end(), cellPointerComp_); -} - -template -template -inline Heap_column::Heap_column( - index columnIndex, const Container_type &nonZeroRowIndices, dimension_type dimension, Row_container_type *rowContainer, Field_operators* operators, Cell_constructor* cellConstructor) - : dim_opt(dimension), - chain_opt([&]{ - if constexpr (Master_matrix::Option_list::is_z2){ - return nonZeroRowIndices.begin() == nonZeroRowIndices.end() ? -1 : *std::prev(nonZeroRowIndices.end()); - } else { - return nonZeroRowIndices.begin() == nonZeroRowIndices.end() ? -1 : std::prev(nonZeroRowIndices.end())->first; - } - }()), - column_(nonZeroRowIndices.size(), nullptr), - insertsSinceLastPrune_(0), - operators_(operators), - cellPool_(cellConstructor) -{ - _verifyCellConstructor(); - - index i = 0; - if constexpr (Master_matrix::Option_list::is_z2){ - for (id_index id : nonZeroRowIndices){ - column_[i++] = cellPool_->construct(id); - } - } else { - for (const auto& p : nonZeroRowIndices){ - column_[i] = cellPool_->construct(p.first); - column_[i++]->set_element(operators_->get_value(p.second)); - } - } - std::make_heap(column_.begin(), column_.end(), cellPointerComp_); -} - -template -template -inline Heap_column::Heap_column( - const Heap_column &column, index columnIndex, Row_container_type *rowContainer, Field_operators* operators, Cell_constructor* cellConstructor) - : dim_opt(static_cast(column)), - chain_opt(static_cast(column)), - column_(column.column_.size(), nullptr), - insertsSinceLastPrune_(0), - operators_(operators == nullptr ? column.operators_ : operators), - cellPool_(cellConstructor == nullptr ? column.cellPool_ : cellConstructor) -{ - index i = 0; - for (const Cell* cell : column.column_){ - if constexpr (Master_matrix::Option_list::is_z2){ - column_[i++] = cellPool_->construct(cell->get_row_index()); - } else { - column_[i] = cellPool_->construct(cell->get_row_index()); - column_[i++]->set_element(cell->get_element()); - } - } - //column.column_ already ordered as a heap, so no need of make_heap. -} - -template -inline Heap_column::~Heap_column() -{ - for (auto* cell : column_){ - cellPool_->destroy(cell); - } -} - -template -inline std::vector::Field_element_type> -Heap_column::get_content(int columnLength) const -{ - bool pivotLength = (columnLength < 0); - if (columnLength < 0 && column_.size() > 0) columnLength = column_.front()->get_row_index() + 1; - else if (columnLength < 0) return std::vector(); - - std::vector container(columnLength, 0); - for (auto it = column_.begin(); it != column_.end(); ++it){ - if ((*it)->get_row_index() < static_cast(columnLength)){ - if constexpr (Master_matrix::Option_list::is_z2){ - container[(*it)->get_row_index()] = !container[(*it)->get_row_index()]; - } else { - container[(*it)->get_row_index()] = operators_->add(container[(*it)->get_row_index()], (*it)->get_element()); - } - } - } - - if (pivotLength){ - while (!container.empty() && container.back() == 0u) container.pop_back(); - } - - return container; -} - -template -inline bool Heap_column::is_non_zero(id_index rowIndex) const -{ - if constexpr (Master_matrix::Option_list::is_z2){ - bool c = false; - for (const Cell* cell : column_){ - if (cell->get_row_index() == rowIndex) c = !c; - } - return c; - } else { - Field_element_type c(0); - for (const Cell* cell : column_){ - if (cell->get_row_index() == rowIndex) c = operators_->add(c, cell->get_element()); - } - return c != Field_operators::get_additive_identity(); - } -} - -template -inline bool Heap_column::is_empty() -{ - Cell* pivot = _pop_pivot(); - if (pivot != nullptr){ - column_.push_back(pivot); - std::push_heap(column_.begin(), column_.end(), cellPointerComp_); - return false; - } - return true; -} - -template -inline std::size_t Heap_column::size() const { - return column_.size(); -} - -template -template -inline void Heap_column::reorder(const Map_type &valueMap, [[maybe_unused]] index columnIndex) -{ - static_assert(!Master_matrix::isNonBasic || Master_matrix::Option_list::is_of_boundary_type, - "Method not available for chain columns."); +template +template +inline Heap_column::Heap_column(index columnIndex, + const Container_type& nonZeroRowIndices, + Row_container_type* rowContainer, + Field_operators* operators, + Cell_constructor* cellConstructor) + : dim_opt(nonZeroRowIndices.size() == 0 ? 0 : nonZeroRowIndices.size() - 1), + chain_opt([&] { + if constexpr (Master_matrix::Option_list::is_z2) { + return nonZeroRowIndices.begin() == nonZeroRowIndices.end() ? -1 : *std::prev(nonZeroRowIndices.end()); + } else { + return nonZeroRowIndices.begin() == nonZeroRowIndices.end() ? -1 : std::prev(nonZeroRowIndices.end())->first; + } + }()), + column_(nonZeroRowIndices.size(), nullptr), + insertsSinceLastPrune_(0), + operators_(operators), + cellPool_(cellConstructor) +{ + static_assert(!Master_matrix::isNonBasic || Master_matrix::Option_list::is_of_boundary_type, + "Constructor not available for chain columns, please specify the dimension of the chain."); + + _verifyCellConstructor(); + + index i = 0; + if constexpr (Master_matrix::Option_list::is_z2) { + for (id_index id : nonZeroRowIndices) { + column_[i++] = cellPool_->construct(id); + } + } else { + for (const auto& p : nonZeroRowIndices) { + column_[i] = cellPool_->construct(p.first); + column_[i++]->set_element(operators_->get_value(p.second)); + } + } + std::make_heap(column_.begin(), column_.end(), cellPointerComp_); +} - Column_type tempCol; - Cell* pivot = _pop_pivot(); - while (pivot != nullptr) { - pivot->set_row_index(valueMap.at(pivot->get_row_index())); - tempCol.push_back(pivot); - pivot = _pop_pivot(); - } - column_.swap(tempCol); - std::make_heap(column_.begin(), column_.end(), cellPointerComp_); +template +template +inline Heap_column::Heap_column( + index columnIndex, + const Container_type& nonZeroRowIndices, + dimension_type dimension, + Row_container_type* rowContainer, + Field_operators* operators, + Cell_constructor* cellConstructor) + : dim_opt(dimension), + chain_opt([&] { + if constexpr (Master_matrix::Option_list::is_z2) { + return nonZeroRowIndices.begin() == nonZeroRowIndices.end() ? -1 : *std::prev(nonZeroRowIndices.end()); + } else { + return nonZeroRowIndices.begin() == nonZeroRowIndices.end() ? -1 : std::prev(nonZeroRowIndices.end())->first; + } + }()), + column_(nonZeroRowIndices.size(), nullptr), + insertsSinceLastPrune_(0), + operators_(operators), + cellPool_(cellConstructor) +{ + _verifyCellConstructor(); + + index i = 0; + if constexpr (Master_matrix::Option_list::is_z2) { + for (id_index id : nonZeroRowIndices) { + column_[i++] = cellPool_->construct(id); + } + } else { + for (const auto& p : nonZeroRowIndices) { + column_[i] = cellPool_->construct(p.first); + column_[i++]->set_element(operators_->get_value(p.second)); + } + } + std::make_heap(column_.begin(), column_.end(), cellPointerComp_); +} - insertsSinceLastPrune_ = 0; +template +template +inline Heap_column::Heap_column(const Heap_column& column, index columnIndex, + Row_container_type* rowContainer, + Field_operators* operators, + Cell_constructor* cellConstructor) + : dim_opt(static_cast(column)), + chain_opt(static_cast(column)), + column_(column.column_.size(), nullptr), + insertsSinceLastPrune_(0), + operators_(operators == nullptr ? column.operators_ : operators), + cellPool_(cellConstructor == nullptr ? column.cellPool_ : cellConstructor) +{ + index i = 0; + for (const Cell* cell : column.column_) { + if constexpr (Master_matrix::Option_list::is_z2) { + column_[i++] = cellPool_->construct(cell->get_row_index()); + } else { + column_[i] = cellPool_->construct(cell->get_row_index()); + column_[i++]->set_element(cell->get_element()); + } + } + // column.column_ already ordered as a heap, so no need of make_heap. } -template -inline void Heap_column::clear() +template +inline Heap_column::~Heap_column() { - static_assert(!Master_matrix::isNonBasic || Master_matrix::Option_list::is_of_boundary_type, - "Method not available for chain columns as a base element should not be empty."); + for (auto* cell : column_) { + cellPool_->destroy(cell); + } +} - for (auto* cell : column_){ - cellPool_->destroy(cell); - } +template +inline std::vector::Field_element_type> +Heap_column::get_content(int columnLength) const +{ + bool pivotLength = (columnLength < 0); + if (columnLength < 0 && column_.size() > 0) + columnLength = column_.front()->get_row_index() + 1; + else if (columnLength < 0) + return std::vector(); + + std::vector container(columnLength, 0); + for (auto it = column_.begin(); it != column_.end(); ++it) { + if ((*it)->get_row_index() < static_cast(columnLength)) { + if constexpr (Master_matrix::Option_list::is_z2) { + container[(*it)->get_row_index()] = !container[(*it)->get_row_index()]; + } else { + container[(*it)->get_row_index()] = operators_->add(container[(*it)->get_row_index()], (*it)->get_element()); + } + } + } + + if (pivotLength) { + while (!container.empty() && container.back() == 0u) container.pop_back(); + } + + return container; +} - column_.clear(); - insertsSinceLastPrune_ = 0; +template +inline bool Heap_column::is_non_zero(id_index rowIndex) const +{ + if constexpr (Master_matrix::Option_list::is_z2) { + bool c = false; + for (const Cell* cell : column_) { + if (cell->get_row_index() == rowIndex) c = !c; + } + return c; + } else { + Field_element_type c(0); + for (const Cell* cell : column_) { + if (cell->get_row_index() == rowIndex) c = operators_->add(c, cell->get_element()); + } + return c != Field_operators::get_additive_identity(); + } } -template -inline void Heap_column::clear(id_index rowIndex) +template +inline bool Heap_column::is_empty() { - static_assert(!Master_matrix::isNonBasic || Master_matrix::Option_list::is_of_boundary_type, - "Method not available for chain columns."); + Cell* pivot = _pop_pivot(); + if (pivot != nullptr) { + column_.push_back(pivot); + std::push_heap(column_.begin(), column_.end(), cellPointerComp_); + return false; + } + return true; +} - column_.erase(std::remove_if(column_.begin(), - column_.end(), - [rowIndex](const Cell* c){ return c->get_row_index() == rowIndex; }), - column_.end()); - std::make_heap(column_.begin(), column_.end(), cellPointerComp_); +template +inline std::size_t Heap_column::size() const +{ + return column_.size(); +} + +template +template +inline void Heap_column::reorder(const Map_type& valueMap, + [[maybe_unused]] index columnIndex) +{ + static_assert(!Master_matrix::isNonBasic || Master_matrix::Option_list::is_of_boundary_type, + "Method not available for chain columns."); + + Column_type tempCol; + Cell* pivot = _pop_pivot(); + while (pivot != nullptr) { + pivot->set_row_index(valueMap.at(pivot->get_row_index())); + tempCol.push_back(pivot); + pivot = _pop_pivot(); + } + column_.swap(tempCol); + std::make_heap(column_.begin(), column_.end(), cellPointerComp_); + + insertsSinceLastPrune_ = 0; } -template -inline typename Heap_column::id_index Heap_column::get_pivot() +template +inline void Heap_column::clear() { - static_assert(Master_matrix::isNonBasic, "Method not available for base columns."); //could technically be, but is the notion usefull then? + static_assert(!Master_matrix::isNonBasic || Master_matrix::Option_list::is_of_boundary_type, + "Method not available for chain columns as a base element should not be empty."); - if constexpr (Master_matrix::Option_list::is_of_boundary_type){ - Cell* pivot = _pop_pivot(); - if (pivot != nullptr){ - column_.push_back(pivot); - std::push_heap(column_.begin(), column_.end(), cellPointerComp_); - return pivot->get_row_index(); - } - return -1; - } else { - return chain_opt::get_pivot(); - } + for (auto* cell : column_) { + cellPool_->destroy(cell); + } + + column_.clear(); + insertsSinceLastPrune_ = 0; } -template -inline typename Heap_column::Field_element_type Heap_column::get_pivot_value() +template +inline void Heap_column::clear(id_index rowIndex) { - static_assert(Master_matrix::isNonBasic, "Method not available for base columns."); //could technically be, but is the notion usefull then? + static_assert(!Master_matrix::isNonBasic || Master_matrix::Option_list::is_of_boundary_type, + "Method not available for chain columns."); - if constexpr (Master_matrix::Option_list::is_z2){ - return 1; - } else { - if constexpr (Master_matrix::Option_list::is_of_boundary_type){ - Cell* pivot = _pop_pivot(); - if (pivot != nullptr){ - column_.push_back(pivot); - std::push_heap(column_.begin(), column_.end(), cellPointerComp_); - return pivot->get_element(); - } - return 0; - } else { - Field_element_type sum(0); - if (chain_opt::get_pivot() == -1) return sum; - for (const Cell* cell : column_){ - if (cell->get_row_index() == chain_opt::get_pivot()) sum = operators_->add(sum, cell->get_element()); - } - return sum; //should not be 0 if properly used. - } - } + column_.erase(std::remove_if(column_.begin(), column_.end(), + [rowIndex](const Cell* c) { return c->get_row_index() == rowIndex; }), + column_.end()); + std::make_heap(column_.begin(), column_.end(), cellPointerComp_); } -template -inline typename Heap_column::iterator -Heap_column::begin() noexcept +template +inline typename Heap_column::id_index +Heap_column::get_pivot() { - return column_.begin(); + static_assert(Master_matrix::isNonBasic, + "Method not available for base columns."); // could technically be, but is the notion usefull then? + + if constexpr (Master_matrix::Option_list::is_of_boundary_type) { + Cell* pivot = _pop_pivot(); + if (pivot != nullptr) { + column_.push_back(pivot); + std::push_heap(column_.begin(), column_.end(), cellPointerComp_); + return pivot->get_row_index(); + } + return -1; + } else { + return chain_opt::get_pivot(); + } } -template -inline typename Heap_column::const_iterator -Heap_column::begin() const noexcept +template +inline typename Heap_column::Field_element_type +Heap_column::get_pivot_value() { - return column_.begin(); + static_assert(Master_matrix::isNonBasic, + "Method not available for base columns."); // could technically be, but is the notion usefull then? + + if constexpr (Master_matrix::Option_list::is_z2) { + return 1; + } else { + if constexpr (Master_matrix::Option_list::is_of_boundary_type) { + Cell* pivot = _pop_pivot(); + if (pivot != nullptr) { + column_.push_back(pivot); + std::push_heap(column_.begin(), column_.end(), cellPointerComp_); + return pivot->get_element(); + } + return 0; + } else { + Field_element_type sum(0); + if (chain_opt::get_pivot() == -1) return sum; + for (const Cell* cell : column_) { + if (cell->get_row_index() == chain_opt::get_pivot()) sum = operators_->add(sum, cell->get_element()); + } + return sum; // should not be 0 if properly used. + } + } } -template -inline typename Heap_column::iterator -Heap_column::end() noexcept +template +inline typename Heap_column::iterator +Heap_column::begin() noexcept { - return column_.end(); + return column_.begin(); } -template -inline typename Heap_column::const_iterator -Heap_column::end() const noexcept +template +inline typename Heap_column::const_iterator +Heap_column::begin() const noexcept { - return column_.end(); + return column_.begin(); } -template -inline typename Heap_column::reverse_iterator -Heap_column::rbegin() noexcept +template +inline typename Heap_column::iterator +Heap_column::end() noexcept { - return column_.rbegin(); + return column_.end(); } -template -inline typename Heap_column::const_reverse_iterator -Heap_column::rbegin() const noexcept +template +inline typename Heap_column::const_iterator +Heap_column::end() const noexcept { - return column_.rbegin(); + return column_.end(); } -template -inline typename Heap_column::reverse_iterator -Heap_column::rend() noexcept +template +inline typename Heap_column::reverse_iterator +Heap_column::rbegin() noexcept { - return column_.rend(); + return column_.rbegin(); } -template -inline typename Heap_column::const_reverse_iterator -Heap_column::rend() const noexcept +template +inline typename Heap_column::const_reverse_iterator +Heap_column::rbegin() const noexcept { - return column_.rend(); + return column_.rbegin(); } -template -template -inline Heap_column & -Heap_column::operator+=(const Cell_range &column) +template +inline typename Heap_column::reverse_iterator +Heap_column::rend() noexcept { - static_assert((!Master_matrix::isNonBasic || std::is_same_v), - "For boundary columns, the range has to be a column of same type to help ensure the validity of the base element."); //could be removed, if we give the responsability to the user. - static_assert((!Master_matrix::isNonBasic || Master_matrix::Option_list::is_of_boundary_type), - "For chain columns, the given column cannot be constant."); + return column_.rend(); +} - _add(column); +template +inline typename Heap_column::const_reverse_iterator +Heap_column::rend() const noexcept +{ + return column_.rend(); +} - return *this; +template +template +inline Heap_column& Heap_column::operator+=( + const Cell_range& column) +{ + static_assert((!Master_matrix::isNonBasic || std::is_same_v), + "For boundary columns, the range has to be a column of same type to help ensure the validity of the " + "base element."); // could be removed, if we give the responsability to the user. + static_assert((!Master_matrix::isNonBasic || Master_matrix::Option_list::is_of_boundary_type), + "For chain columns, the given column cannot be constant."); + + _add(column); + + return *this; +} + +template +inline Heap_column& Heap_column::operator+=( + Heap_column& column) +{ + if constexpr (Master_matrix::isNonBasic && !Master_matrix::Option_list::is_of_boundary_type) { + // assumes that the addition never zeros out this column. + if (_add(column)) { + chain_opt::swap_pivots(column); + dim_opt::swap_dimension(column); + } + } else { + _add(column); + } + + return *this; +} + +template +inline Heap_column& Heap_column::operator*=( + unsigned int v) +{ + if constexpr (Master_matrix::Option_list::is_z2) { + if (v % 2 == 0) { + if constexpr (Master_matrix::isNonBasic && !Master_matrix::Option_list::is_of_boundary_type) { + throw std::invalid_argument("A chain column should not be multiplied by 0."); + } else { + clear(); + } + } + } else { + Field_element_type val = operators_->get_value(v); + + if (val == Field_operators::get_additive_identity()) { + if constexpr (Master_matrix::isNonBasic && !Master_matrix::Option_list::is_of_boundary_type) { + throw std::invalid_argument("A chain column should not be multiplied by 0."); + } else { + clear(); + } + return *this; + } + + if (val == Field_operators::get_multiplicative_identity()) return *this; + + for (Cell* cell : column_) { + cell->get_element() = operators_->multiply(cell->get_element(), val); + } + } + + return *this; +} + +template +template +inline Heap_column& Heap_column::multiply_and_add( + const Field_element_type& val, const Cell_range& column) +{ + static_assert((!Master_matrix::isNonBasic || std::is_same_v), + "For boundary columns, the range has to be a column of same type to help ensure the validity of the " + "base element."); // could be removed, if we give the responsability to the user. + static_assert((!Master_matrix::isNonBasic || Master_matrix::Option_list::is_of_boundary_type), + "For chain columns, the given column cannot be constant."); + + if constexpr (Master_matrix::Option_list::is_z2) { + if (val) { + _add(column); + } else { + clear(); + _add(column); + } + } else { + _multiply_and_add(val, column); + } + + return *this; +} + +template +inline Heap_column& Heap_column::multiply_and_add( + const Field_element_type& val, Heap_column& column) +{ + if constexpr (Master_matrix::isNonBasic && !Master_matrix::Option_list::is_of_boundary_type) { + // assumes that the addition never zeros out this column. + if constexpr (Master_matrix::Option_list::is_z2) { + if (val) { + if (_add(column)) { + chain_opt::swap_pivots(column); + dim_opt::swap_dimension(column); + } + } else { + throw std::invalid_argument("A chain column should not be multiplied by 0."); + } + } else { + if (_multiply_and_add(val, column)) { + chain_opt::swap_pivots(column); + dim_opt::swap_dimension(column); + } + } + } else { + if constexpr (Master_matrix::Option_list::is_z2) { + if (val) { + _add(column); + } else { + clear(); + _add(column); + } + } else { + _multiply_and_add(val, column); + } + } + + return *this; +} + +template +template +inline Heap_column& Heap_column::multiply_and_add( + const Cell_range& column, const Field_element_type& val) +{ + static_assert((!Master_matrix::isNonBasic || std::is_same_v), + "For boundary columns, the range has to be a column of same type to help ensure the validity of the " + "base element."); // could be removed, if we give the responsability to the user. + static_assert((!Master_matrix::isNonBasic || Master_matrix::Option_list::is_of_boundary_type), + "For chain columns, the given column cannot be constant."); + + if constexpr (Master_matrix::Option_list::is_z2) { + if (val) { + _add(column); + } + } else { + _multiply_and_add(column, val); + } + + return *this; } -template -inline Heap_column & -Heap_column::operator+=(Heap_column &column) +template +inline Heap_column& Heap_column::multiply_and_add( + Heap_column& column, const Field_element_type& val) { - if constexpr (Master_matrix::isNonBasic && !Master_matrix::Option_list::is_of_boundary_type){ - //assumes that the addition never zeros out this column. - if (_add(column)){ - chain_opt::swap_pivots(column); - dim_opt::swap_dimension(column); - } - } else { - _add(column); - } + if constexpr (Master_matrix::isNonBasic && !Master_matrix::Option_list::is_of_boundary_type) { + // assumes that the addition never zeros out this column. + if constexpr (Master_matrix::Option_list::is_z2) { + if (val) { + if (_add(column)) { + chain_opt::swap_pivots(column); + dim_opt::swap_dimension(column); + } + } + } else { + if (_multiply_and_add(column, val)) { + chain_opt::swap_pivots(column); + dim_opt::swap_dimension(column); + } + } + } else { + if constexpr (Master_matrix::Option_list::is_z2) { + if (val) { + _add(column); + } + } else { + _multiply_and_add(column, val); + } + } + + return *this; +} - return *this; +template +inline Heap_column& Heap_column::operator=( + const Heap_column& other) +{ + static_assert(!Master_matrix::Option_list::has_row_access, "= assignement not enabled with row access option."); + + dim_opt::operator=(other); + chain_opt::operator=(other); + + while (column_.size() > other.column_.size()) { + if (column_.back() != nullptr) cellPool_->destroy(column_.back()); + column_.pop_back(); + } + + column_.resize(other.column_.size(), nullptr); + index i = 0; + for (const Cell* cell : other.column_) { + if (column_[i] != nullptr) { + cellPool_->destroy(column_[i]); + } + if constexpr (Master_matrix::Option_list::is_z2) { + column_[i++] = other.cellPool_->construct(cell->get_row_index()); + } else { + column_[i] = other.cellPool_->construct(cell->get_row_index()); + column_[i++]->set_element(cell->get_element()); + } + } + insertsSinceLastPrune_ = other.insertsSinceLastPrune_; + operators_ = other.operators_; + cellPool_ = other.cellPool_; + + return *this; } -template -inline Heap_column & -Heap_column::operator*=(unsigned int v) +template +inline void Heap_column::_prune() { - if constexpr (Master_matrix::Option_list::is_z2){ - if (v % 2 == 0){ - if constexpr (Master_matrix::isNonBasic && !Master_matrix::Option_list::is_of_boundary_type){ - throw std::invalid_argument("A chain column should not be multiplied by 0."); - } else { - clear(); - } - } - } else { - Field_element_type val = operators_->get_value(v); + if (insertsSinceLastPrune_ == 0) return; + + Column_type tempCol; + Cell* pivot = _pop_pivot(); + while (pivot != nullptr) { + tempCol.push_back(pivot); + pivot = _pop_pivot(); + } + column_.swap(tempCol); + std::make_heap(column_.begin(), column_.end(), cellPointerComp_); + + insertsSinceLastPrune_ = 0; +} - if (val == Field_operators::get_additive_identity()) { - if constexpr (Master_matrix::isNonBasic && !Master_matrix::Option_list::is_of_boundary_type){ - throw std::invalid_argument("A chain column should not be multiplied by 0."); - } else { - clear(); - } - return *this; - } +template +inline typename Heap_column::Cell* +Heap_column::_pop_pivot() +{ + if (column_.empty()) { + return nullptr; + } + + Cell* pivot = column_.front(); + std::pop_heap(column_.begin(), column_.end(), cellPointerComp_); + column_.pop_back(); + if constexpr (Master_matrix::Option_list::is_z2) { + while (!column_.empty() && column_.front()->get_row_index() == pivot->get_row_index()) { + std::pop_heap(column_.begin(), column_.end(), cellPointerComp_); + column_.pop_back(); + + if (column_.empty()) { + return nullptr; + } + pivot = column_.front(); + std::pop_heap(column_.begin(), column_.end(), cellPointerComp_); + column_.pop_back(); + } + } else { + while (!column_.empty() && column_.front()->get_row_index() == pivot->get_row_index()) { + pivot->get_element() = operators_->add(pivot->get_element(), column_.front()->get_element()); + std::pop_heap(column_.begin(), column_.end(), cellPointerComp_); + column_.pop_back(); + } + + if (pivot->get_element() == Field_operators::get_additive_identity()) return _pop_pivot(); + } + + return pivot; +} - if (val == Field_operators::get_multiplicative_identity()) return *this; +template +template +inline bool Heap_column::_add(const Cell_range& column) +{ + if (column.begin() == column.end()) return false; + if (column_.empty()) { // chain should never enter here. + column_.resize(column.size()); + index i = 0; + for (const Cell& cell : column) { + if constexpr (Master_matrix::Option_list::is_z2) { + column_[i++] = cellPool_->construct(cell.get_row_index()); + } else { + column_[i] = cellPool_->construct(cell.get_row_index()); + column_[i++]->set_element(cell.get_element()); + } + } + insertsSinceLastPrune_ = column_.size(); + return true; + } + + Field_element_type pivotVal(1); + + if constexpr (Master_matrix::isNonBasic && !Master_matrix::Option_list::is_of_boundary_type) + pivotVal = get_pivot_value(); + + for (const Cell& cell : column) { + ++insertsSinceLastPrune_; + if constexpr (Master_matrix::Option_list::is_z2) { + if constexpr (Master_matrix::isNonBasic && !Master_matrix::Option_list::is_of_boundary_type) { + if (cell.get_row_index() == chain_opt::get_pivot()) pivotVal = !pivotVal; + } + column_.push_back(cellPool_->construct(cell.get_row_index())); + } else { + if constexpr (Master_matrix::isNonBasic && !Master_matrix::Option_list::is_of_boundary_type) { + if (cell.get_row_index() == chain_opt::get_pivot()) pivotVal = operators_->add(pivotVal, cell.get_element()); + } + column_.push_back(cellPool_->construct(cell.get_row_index())); + column_.back()->set_element(cell.get_element()); + } + std::push_heap(column_.begin(), column_.end(), cellPointerComp_); + } + + if (2 * insertsSinceLastPrune_ > column_.size()) _prune(); + + return pivotVal == Field_operators::get_additive_identity(); +} - for (Cell* cell : column_){ - cell->get_element() = operators_->multiply(cell->get_element(), val); - } - } - - return *this; -} - -template -template -inline Heap_column & -Heap_column::multiply_and_add(const Field_element_type& val, const Cell_range& column) -{ - static_assert((!Master_matrix::isNonBasic || std::is_same_v), - "For boundary columns, the range has to be a column of same type to help ensure the validity of the base element."); //could be removed, if we give the responsability to the user. - static_assert((!Master_matrix::isNonBasic || Master_matrix::Option_list::is_of_boundary_type), - "For chain columns, the given column cannot be constant."); - - if constexpr (Master_matrix::Option_list::is_z2){ - if (val){ - _add(column); - } else { - clear(); - _add(column); - } - } else { - _multiply_and_add(val, column); - } - - return *this; -} - -template -inline Heap_column & -Heap_column::multiply_and_add(const Field_element_type& val, Heap_column& column) -{ - if constexpr (Master_matrix::isNonBasic && !Master_matrix::Option_list::is_of_boundary_type){ - //assumes that the addition never zeros out this column. - if constexpr (Master_matrix::Option_list::is_z2){ - if (val){ - if (_add(column)) { - chain_opt::swap_pivots(column); - dim_opt::swap_dimension(column); - } - } else { - throw std::invalid_argument("A chain column should not be multiplied by 0."); - } - } else { - if (_multiply_and_add(val, column)) { - chain_opt::swap_pivots(column); - dim_opt::swap_dimension(column); - } - } - } else { - if constexpr (Master_matrix::Option_list::is_z2){ - if (val){ - _add(column); - } else { - clear(); - _add(column); - } - } else { - _multiply_and_add(val, column); - } - } - - return *this; -} - -template -template -inline Heap_column & -Heap_column::multiply_and_add(const Cell_range& column, const Field_element_type& val) -{ - static_assert((!Master_matrix::isNonBasic || std::is_same_v), - "For boundary columns, the range has to be a column of same type to help ensure the validity of the base element."); //could be removed, if we give the responsability to the user. - static_assert((!Master_matrix::isNonBasic || Master_matrix::Option_list::is_of_boundary_type), - "For chain columns, the given column cannot be constant."); - - if constexpr (Master_matrix::Option_list::is_z2){ - if (val){ - _add(column); - } - } else { - _multiply_and_add(column, val); - } - - return *this; -} - -template -inline Heap_column & -Heap_column::multiply_and_add(Heap_column& column, const Field_element_type& val) -{ - if constexpr (Master_matrix::isNonBasic && !Master_matrix::Option_list::is_of_boundary_type){ - //assumes that the addition never zeros out this column. - if constexpr (Master_matrix::Option_list::is_z2){ - if (val){ - if (_add(column)){ - chain_opt::swap_pivots(column); - dim_opt::swap_dimension(column); - } - } - } else { - if (_multiply_and_add(column, val)){ - chain_opt::swap_pivots(column); - dim_opt::swap_dimension(column); - } - } - } else { - if constexpr (Master_matrix::Option_list::is_z2){ - if (val){ - _add(column); - } - } else { - _multiply_and_add(column, val); - } - } - - return *this; -} - -template -inline Heap_column & -Heap_column::operator=(const Heap_column& other) -{ - static_assert(!Master_matrix::Option_list::has_row_access, "= assignement not enabled with row access option."); - - dim_opt::operator=(other); - chain_opt::operator=(other); - - while (column_.size() > other.column_.size()) { - if (column_.back() != nullptr) cellPool_->destroy(column_.back()); - column_.pop_back(); - } - - column_.resize(other.column_.size(), nullptr); - index i = 0; - for (const Cell* cell : other.column_){ - if (column_[i] != nullptr){ - cellPool_->destroy(column_[i]); - } - if constexpr (Master_matrix::Option_list::is_z2){ - column_[i++] = other.cellPool_->construct(cell->get_row_index()); - } else { - column_[i] = other.cellPool_->construct(cell->get_row_index()); - column_[i++]->set_element(cell->get_element()); - } - } - insertsSinceLastPrune_ = other.insertsSinceLastPrune_; - operators_ = other.operators_; - cellPool_ = other.cellPool_; - - return *this; -} - -template -inline void Heap_column::_prune() -{ - if (insertsSinceLastPrune_ == 0) return; - - Column_type tempCol; - Cell* pivot = _pop_pivot(); - while (pivot != nullptr) { - tempCol.push_back(pivot); - pivot = _pop_pivot(); - } - column_.swap(tempCol); - std::make_heap(column_.begin(), column_.end(), cellPointerComp_); - - insertsSinceLastPrune_ = 0; -} - -template -inline typename Heap_column::Cell* Heap_column::_pop_pivot() -{ - if (column_.empty()) { - return nullptr; - } - - Cell* pivot = column_.front(); - std::pop_heap(column_.begin(), column_.end(), cellPointerComp_); - column_.pop_back(); - if constexpr (Master_matrix::Option_list::is_z2){ - while (!column_.empty() && column_.front()->get_row_index() == pivot->get_row_index()) - { - std::pop_heap(column_.begin(), column_.end(), cellPointerComp_); - column_.pop_back(); - - if (column_.empty()) { - return nullptr; - } - pivot = column_.front(); - std::pop_heap(column_.begin(), column_.end(), cellPointerComp_); - column_.pop_back(); - } - } else { - while (!column_.empty() && column_.front()->get_row_index() == pivot->get_row_index()) - { - pivot->get_element() = operators_->add(pivot->get_element(), column_.front()->get_element()); - std::pop_heap(column_.begin(), column_.end(), cellPointerComp_); - column_.pop_back(); - } - - if (pivot->get_element() == Field_operators::get_additive_identity()) return _pop_pivot(); - } - - return pivot; -} - -template -template -inline bool Heap_column::_add(const Cell_range &column) -{ - if (column.begin() == column.end()) return false; - if (column_.empty()){ //chain should never enter here. - column_.resize(column.size()); - index i = 0; - for (const Cell& cell : column){ - if constexpr (Master_matrix::Option_list::is_z2){ - column_[i++] = cellPool_->construct(cell.get_row_index()); - } else { - column_[i] = cellPool_->construct(cell.get_row_index()); - column_[i++]->set_element(cell.get_element()); - } - } - insertsSinceLastPrune_ = column_.size(); - return true; - } - - Field_element_type pivotVal(1); - - if constexpr (Master_matrix::isNonBasic && !Master_matrix::Option_list::is_of_boundary_type) - pivotVal = get_pivot_value(); - - for (const Cell& cell : column) { - ++insertsSinceLastPrune_; - if constexpr (Master_matrix::Option_list::is_z2){ - if constexpr (Master_matrix::isNonBasic && !Master_matrix::Option_list::is_of_boundary_type){ - if (cell.get_row_index() == chain_opt::get_pivot()) pivotVal = !pivotVal; - } - column_.push_back(cellPool_->construct(cell.get_row_index())); - } else { - if constexpr (Master_matrix::isNonBasic && !Master_matrix::Option_list::is_of_boundary_type){ - if (cell.get_row_index() == chain_opt::get_pivot()) pivotVal = operators_->add(pivotVal, cell.get_element()); - } - column_.push_back(cellPool_->construct(cell.get_row_index())); - column_.back()->set_element(cell.get_element()); - } - std::push_heap(column_.begin(), column_.end(), cellPointerComp_); - } - - if (2 * insertsSinceLastPrune_ > column_.size()) _prune(); - - return pivotVal == Field_operators::get_additive_identity(); -} - -template -template -inline bool Heap_column::_multiply_and_add(const Field_element_type& val, const Cell_range& column) -{ - if (val == 0u) { - if constexpr (Master_matrix::isNonBasic && !Master_matrix::Option_list::is_of_boundary_type){ - throw std::invalid_argument("A chain column should not be multiplied by 0."); - //this would not only mess up the base, but also the pivots stored. - } else { - clear(); - } - } - if (column_.empty()){ //chain should never enter here. - column_.resize(column.size()); - index i = 0; - for (const Cell& cell : column){ - column_[i] = cellPool_->construct(cell.get_row_index()); - column_[i++]->set_element(cell.get_element()); - } - insertsSinceLastPrune_ = column_.size(); - return true; - } - - Field_element_type pivotVal(0); - - for (Cell* cell : column_){ - cell->get_element() = operators_->multiply(cell->get_element(), val); - if constexpr (Master_matrix::isNonBasic && !Master_matrix::Option_list::is_of_boundary_type){ - if (cell->get_row_index() == chain_opt::get_pivot()) pivotVal = operators_->add(pivotVal, cell->get_element()); - } - } - - for (const Cell& cell : column) { - ++insertsSinceLastPrune_; - if constexpr (Master_matrix::isNonBasic && !Master_matrix::Option_list::is_of_boundary_type){ - if (cell.get_row_index() == chain_opt::get_pivot()) pivotVal = operators_->add(pivotVal, cell.get_element()); - } - column_.push_back(cellPool_->construct(cell.get_row_index())); - column_.back()->set_element(cell.get_element()); - std::push_heap(column_.begin(), column_.end(), cellPointerComp_); - } - - if (2 * insertsSinceLastPrune_ > column_.size()) _prune(); - - return pivotVal == Field_operators::get_additive_identity(); -} - -template -template -inline bool Heap_column::_multiply_and_add(const Cell_range& column, const Field_element_type& val) -{ - if (val == 0u || column.begin() == column.end()) { - return false; - } - if (column_.empty()){ //chain should never enter here. - column_.resize(column.size()); - index i = 0; - for (const Cell& cell : column){ - column_[i] = cellPool_->construct(cell.get_row_index()); - column_[i++]->set_element(cell.get_element()); - } - insertsSinceLastPrune_ = column_.size(); - return true; - } - - Field_element_type pivotVal(1); - - if constexpr (Master_matrix::isNonBasic && !Master_matrix::Option_list::is_of_boundary_type) - pivotVal = get_pivot_value(); - - for (const Cell& cell : column) { - ++insertsSinceLastPrune_; - if constexpr (Master_matrix::isNonBasic && !Master_matrix::Option_list::is_of_boundary_type){ - if (cell.get_row_index() == chain_opt::get_pivot()) pivotVal = operators_->multiply_and_add(cell.get_element(), val, pivotVal); - } - column_.push_back(cellPool_->construct(cell.get_row_index())); - column_.back()->set_element(operators_->multiply(cell.get_element(), val)); - std::push_heap(column_.begin(), column_.end(), cellPointerComp_); - } - - if (2 * insertsSinceLastPrune_ > column_.size()) _prune(); - - return pivotVal == 0u; -} - -} //namespace persistence_matrix -} //namespace Gudhi - -template -struct std::hash > -{ - size_t operator()(const Gudhi::persistence_matrix::Heap_column& column) const - { - std::size_t seed = 0; - unsigned int i = 0; - for (bool val : column.get_content()){ - seed ^= std::hash()(i++ * val) + 0x9e3779b9 + (seed << 6) + (seed >> 2); - } - return seed; - } +template +template +inline bool Heap_column::_multiply_and_add(const Field_element_type& val, + const Cell_range& column) +{ + if (val == 0u) { + if constexpr (Master_matrix::isNonBasic && !Master_matrix::Option_list::is_of_boundary_type) { + throw std::invalid_argument("A chain column should not be multiplied by 0."); + // this would not only mess up the base, but also the pivots stored. + } else { + clear(); + } + } + if (column_.empty()) { // chain should never enter here. + column_.resize(column.size()); + index i = 0; + for (const Cell& cell : column) { + column_[i] = cellPool_->construct(cell.get_row_index()); + column_[i++]->set_element(cell.get_element()); + } + insertsSinceLastPrune_ = column_.size(); + return true; + } + + Field_element_type pivotVal(0); + + for (Cell* cell : column_) { + cell->get_element() = operators_->multiply(cell->get_element(), val); + if constexpr (Master_matrix::isNonBasic && !Master_matrix::Option_list::is_of_boundary_type) { + if (cell->get_row_index() == chain_opt::get_pivot()) pivotVal = operators_->add(pivotVal, cell->get_element()); + } + } + + for (const Cell& cell : column) { + ++insertsSinceLastPrune_; + if constexpr (Master_matrix::isNonBasic && !Master_matrix::Option_list::is_of_boundary_type) { + if (cell.get_row_index() == chain_opt::get_pivot()) pivotVal = operators_->add(pivotVal, cell.get_element()); + } + column_.push_back(cellPool_->construct(cell.get_row_index())); + column_.back()->set_element(cell.get_element()); + std::push_heap(column_.begin(), column_.end(), cellPointerComp_); + } + + if (2 * insertsSinceLastPrune_ > column_.size()) _prune(); + + return pivotVal == Field_operators::get_additive_identity(); +} + +template +template +inline bool Heap_column::_multiply_and_add(const Cell_range& column, + const Field_element_type& val) +{ + if (val == 0u || column.begin() == column.end()) { + return false; + } + if (column_.empty()) { // chain should never enter here. + column_.resize(column.size()); + index i = 0; + for (const Cell& cell : column) { + column_[i] = cellPool_->construct(cell.get_row_index()); + column_[i++]->set_element(cell.get_element()); + } + insertsSinceLastPrune_ = column_.size(); + return true; + } + + Field_element_type pivotVal(1); + + if constexpr (Master_matrix::isNonBasic && !Master_matrix::Option_list::is_of_boundary_type) + pivotVal = get_pivot_value(); + + for (const Cell& cell : column) { + ++insertsSinceLastPrune_; + if constexpr (Master_matrix::isNonBasic && !Master_matrix::Option_list::is_of_boundary_type) { + if (cell.get_row_index() == chain_opt::get_pivot()) + pivotVal = operators_->multiply_and_add(cell.get_element(), val, pivotVal); + } + column_.push_back(cellPool_->construct(cell.get_row_index())); + column_.back()->set_element(operators_->multiply(cell.get_element(), val)); + std::push_heap(column_.begin(), column_.end(), cellPointerComp_); + } + + if (2 * insertsSinceLastPrune_ > column_.size()) _prune(); + + return pivotVal == 0u; +} + +} // namespace persistence_matrix +} // namespace Gudhi + +/** + * @brief Hash method for @ref Gudhi::persistence_matrix::Heap_column. + * + * @tparam Master_matrix Template parameter of @ref Gudhi::persistence_matrix::Heap_column. + * @tparam Cell_constructor Template parameter of @ref Gudhi::persistence_matrix::Heap_column. + */ +template +struct std::hash > +{ + size_t operator()(const Gudhi::persistence_matrix::Heap_column& column) const { + std::size_t seed = 0; + unsigned int i = 0; + for (bool val : column.get_content()) { + seed ^= std::hash()(i++ * val) + 0x9e3779b9 + (seed << 6) + (seed >> 2); + } + return seed; + } }; -#endif // PM_HEAP_COLUMN_H +#endif // PM_HEAP_COLUMN_H diff --git a/src/Persistence_matrix/include/gudhi/Persistence_matrix/columns/intrusive_list_column.h b/src/Persistence_matrix/include/gudhi/Persistence_matrix/columns/intrusive_list_column.h index 190781e337..a8691a2055 100644 --- a/src/Persistence_matrix/include/gudhi/Persistence_matrix/columns/intrusive_list_column.h +++ b/src/Persistence_matrix/include/gudhi/Persistence_matrix/columns/intrusive_list_column.h @@ -2,19 +2,26 @@ * 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-24 Inria * * Modification(s): * - YYYY/MM Author: Description of the modification */ +/** + * @file intrusive_list_column.h + * @author Hannah Schreiber + * @brief Contains the @ref Intrusive_list_column class. + * Also defines the std::hash method for @ref Intrusive_list_column. + */ + #ifndef PM_INTRUSIVE_LIST_COLUMN_H #define PM_INTRUSIVE_LIST_COLUMN_H #include #include #include -#include //std::swap, std::move & std::exchange +#include //std::swap, std::move & std::exchange #include @@ -23,965 +30,1012 @@ namespace Gudhi { namespace persistence_matrix { -template > -class Intrusive_list_column : public Master_matrix::Row_access_option, - public Master_matrix::Column_dimension_option, - public Master_matrix::Chain_column_option -{ -public: - using Master = Master_matrix; - using Field_operators = typename Master_matrix::Field_operators; - using Field_element_type = typename Master_matrix::element_type; - using index = typename Master_matrix::index; - using id_index = typename Master_matrix::id_index; - using dimension_type = typename Master_matrix::dimension_type; - - using Cell = typename Master_matrix::Cell_type; - using Column_type = boost::intrusive::list < - Cell, - boost::intrusive::constant_time_size, - boost::intrusive::base_hook< typename Master_matrix::base_hook_matrix_list_column > - >; - using iterator = typename Column_type::iterator; - using const_iterator = typename Column_type::const_iterator; - using reverse_iterator = typename Column_type::reverse_iterator; - using const_reverse_iterator = typename Column_type::const_reverse_iterator; - - Intrusive_list_column(Field_operators* operators = nullptr, Cell_constructor* cellConstructor = nullptr); - template - Intrusive_list_column(const Container_type& nonZeroRowIndices, Field_operators* operators, Cell_constructor* cellConstructor); //has to be a boundary for boundary, has no sense for chain if dimension is needed - template - Intrusive_list_column(index columnIndex, const Container_type& nonZeroRowIndices, Row_container_type *rowContainer, Field_operators* operators, Cell_constructor* cellConstructor); //has to be a boundary for boundary, has no sense for chain if dimension is needed - template - Intrusive_list_column(const Container_type& nonZeroChainRowIndices, dimension_type dimension, Field_operators* operators, Cell_constructor* cellConstructor); //dimension gets ignored for base - template - Intrusive_list_column(index columnIndex, const Container_type& nonZeroChainRowIndices, dimension_type dimension, Row_container_type *rowContainer, Field_operators* operators, Cell_constructor* cellConstructor); //dimension gets ignored for base - Intrusive_list_column(const Intrusive_list_column& column, Field_operators* operators = nullptr, Cell_constructor* cellConstructor = nullptr); - template - Intrusive_list_column(const Intrusive_list_column& column, index columnIndex, Row_container_type *rowContainer, Field_operators* operators = nullptr, Cell_constructor* cellConstructor = nullptr); - Intrusive_list_column(Intrusive_list_column&& column) noexcept; - ~Intrusive_list_column(); - - std::vector get_content(int columnLength = -1) const; - bool is_non_zero(id_index rowIndex) const; - bool is_empty() const; - std::size_t size() const; - - //**************** - //only for base and boundary - template - void reorder(const Map_type& valueMap, [[maybe_unused]] index columnIndex = -1); //used for lazy row swaps - void clear(); - void clear(id_index rowIndex); - //**************** - - //**************** - //only for chain and boundary - id_index get_pivot() const; - Field_element_type get_pivot_value() const; - //**************** - - iterator begin() noexcept; - const_iterator begin() const noexcept; - iterator end() noexcept; - const_iterator end() const noexcept; - reverse_iterator rbegin() noexcept; - const_reverse_iterator rbegin() const noexcept; - reverse_iterator rend() noexcept; - const_reverse_iterator rend() const noexcept; - - template - Intrusive_list_column& operator+=(const Cell_range& column); //for base & boundary except vector - Intrusive_list_column& operator+=(Intrusive_list_column &column); //for chain and vector - - Intrusive_list_column& operator*=(const Field_element_type& val); - - //this = v * this + column - template - Intrusive_list_column& multiply_and_add(const Field_element_type& val, const Cell_range& column); //for base & boundary except vector - Intrusive_list_column& multiply_and_add(const Field_element_type& val, Intrusive_list_column& column); //for chain and vector - //this = this + column * v - template - Intrusive_list_column& multiply_and_add(const Cell_range& column, const Field_element_type& val); //for base & boundary except vector - Intrusive_list_column& multiply_and_add(Intrusive_list_column& column, const Field_element_type& val); //for chain and vector - - friend bool operator==(const Intrusive_list_column& c1, const Intrusive_list_column& c2){ - if (&c1 == &c2) return true; - - if constexpr (Master_matrix::Option_list::is_z2){ - return c1.column_ == c2.column_; - } else { - auto it1 = c1.column_.begin(); - auto it2 = c2.column_.begin(); - if (c1.column_.size() != c2.column_.size()) return false; - while (it1 != c1.column_.end() && it2 != c2.column_.end()) { - if (it1->get_row_index() != it2->get_row_index() || it1->get_element() != it2->get_element()) - return false; - ++it1; ++it2; - } - return true; - } - } - friend bool operator<(const Intrusive_list_column& c1, const Intrusive_list_column& c2){ - if (&c1 == &c2) return false; - - if constexpr (Master_matrix::Option_list::is_z2){ - return c1.column_ < c2.column_; - } else { - auto it1 = c1.column_.begin(); - auto it2 = c2.column_.begin(); - while (it1 != c1.column_.end() && it2 != c2.column_.end()) { - if (it1->get_row_index() != it2->get_row_index()) - return it1->get_row_index() < it2->get_row_index(); - if (it1->get_element() != it2->get_element()) - return it1->get_element() < it2->get_element(); - ++it1; ++it2; - } - return it2 != c2.column_.end(); - } - } - - // void set_operators(Field_operators* operators){ operators_ = operators; } - - //Disabled with row access. - Intrusive_list_column& operator=(const Intrusive_list_column& other); - - friend void swap(Intrusive_list_column& col1, Intrusive_list_column& col2){ - swap(static_cast(col1), - static_cast(col2)); - swap(static_cast(col1), - static_cast(col2)); - swap(static_cast(col1), - static_cast(col2)); - col1.column_.swap(col2.column_); - std::swap(col1.operators_, col2.operators_); - std::swap(col1.cellPool_, col2.cellPool_); - } - -private: - using ra_opt = typename Master_matrix::Row_access_option; - using dim_opt = typename Master_matrix::Column_dimension_option; - using chain_opt = typename Master_matrix::Chain_column_option; - - //Cloner object function for boost intrusive container - struct new_cloner - { - new_cloner(Cell_constructor* cellPool) : cellPool_(cellPool) - {}; - - Cell *operator()(const Cell& clone_this) - { - return cellPool_->construct(clone_this); - } - - Cell_constructor* cellPool_; - }; - - //The disposer object function for boost intrusive container - struct delete_disposer - { - delete_disposer(){}; - delete_disposer(Intrusive_list_column* col) : col_(col) - {}; - - void operator()(Cell *delete_this) - { - if constexpr (Master_matrix::Option_list::has_row_access) - col_->unlink(delete_this); - col_->cellPool_->destroy(delete_this); - } - - Intrusive_list_column* col_; - }; - - Field_operators* operators_; - Cell_constructor* cellPool_; - Column_type column_; - - void _delete_cell(iterator& it); - void _insert_cell(const Field_element_type& value, id_index rowIndex, const iterator& position); - void _insert_cell(id_index rowIndex, const iterator& position); - template - bool _add(const Cell_range& column); - template - bool _multiply_and_add(const Field_element_type& val, const Cell_range& column); - template - bool _multiply_and_add(const Cell_range& column, const Field_element_type& val); - - void _verifyCellConstructor(){ - if (cellPool_ == nullptr){ - if constexpr (std::is_same_v >){ - cellPool_ = &Master_matrix::defaultCellConstructor; - } else { - throw std::invalid_argument("Cell constructor pointer cannot be null."); - } - } - } +/** + * @brief Column class following the @ref PersistenceMatrixColumn concept. + * + * Column based on a intrusive list structure. The cells are always ordered by row index and only non-zero values + * are stored uniquely in the underlying container. + * + * @tparam Master_matrix An instanciation of @ref Matrix from which all types and options are deduced. + * @tparam Cell_constructor Factory of @ref Cell classes. + */ +template > +class Intrusive_list_column : public Master_matrix::Row_access_option, + public Master_matrix::Column_dimension_option, + public Master_matrix::Chain_column_option +{ + public: + using Master = Master_matrix; + using Field_operators = typename Master_matrix::Field_operators; + using Field_element_type = typename Master_matrix::element_type; + using index = typename Master_matrix::index; + using id_index = typename Master_matrix::id_index; + using dimension_type = typename Master_matrix::dimension_type; + + using Cell = typename Master_matrix::Cell_type; + using Column_type = + boost::intrusive::list, + boost::intrusive::base_hook >; + using iterator = typename Column_type::iterator; + using const_iterator = typename Column_type::const_iterator; + using reverse_iterator = typename Column_type::reverse_iterator; + using const_reverse_iterator = typename Column_type::const_reverse_iterator; + + Intrusive_list_column(Field_operators* operators = nullptr, Cell_constructor* cellConstructor = nullptr); + template + Intrusive_list_column(const Container_type& nonZeroRowIndices, + Field_operators* operators, + Cell_constructor* cellConstructor); + template + Intrusive_list_column(index columnIndex, + const Container_type& nonZeroRowIndices, + Row_container_type* rowContainer, + Field_operators* operators, + Cell_constructor* cellConstructor); + template + Intrusive_list_column(const Container_type& nonZeroChainRowIndices, + dimension_type dimension, + Field_operators* operators, + Cell_constructor* cellConstructor); + template + Intrusive_list_column(index columnIndex, + const Container_type& nonZeroChainRowIndices, + dimension_type dimension, + Row_container_type* rowContainer, + Field_operators* operators, + Cell_constructor* cellConstructor); + Intrusive_list_column(const Intrusive_list_column& column, + Field_operators* operators = nullptr, + Cell_constructor* cellConstructor = nullptr); + template + Intrusive_list_column(const Intrusive_list_column& column, + index columnIndex, + Row_container_type* rowContainer, + Field_operators* operators = nullptr, + Cell_constructor* cellConstructor = nullptr); + Intrusive_list_column(Intrusive_list_column&& column) noexcept; + ~Intrusive_list_column(); + + std::vector get_content(int columnLength = -1) const; + bool is_non_zero(id_index rowIndex) const; + bool is_empty() const; + std::size_t size() const; + + template + void reorder(const Map_type& valueMap, [[maybe_unused]] index columnIndex = -1); + void clear(); + void clear(id_index rowIndex); + + id_index get_pivot() const; + Field_element_type get_pivot_value() const; + + iterator begin() noexcept; + const_iterator begin() const noexcept; + iterator end() noexcept; + const_iterator end() const noexcept; + reverse_iterator rbegin() noexcept; + const_reverse_iterator rbegin() const noexcept; + reverse_iterator rend() noexcept; + const_reverse_iterator rend() const noexcept; + + template + Intrusive_list_column& operator+=(const Cell_range& column); + Intrusive_list_column& operator+=(Intrusive_list_column& column); + + Intrusive_list_column& operator*=(const Field_element_type& val); + + // this = v * this + column + template + Intrusive_list_column& multiply_and_add(const Field_element_type& val, const Cell_range& column); + Intrusive_list_column& multiply_and_add(const Field_element_type& val, Intrusive_list_column& column); + // this = this + column * v + template + Intrusive_list_column& multiply_and_add(const Cell_range& column, const Field_element_type& val); + Intrusive_list_column& multiply_and_add(Intrusive_list_column& column, const Field_element_type& val); + + friend bool operator==(const Intrusive_list_column& c1, const Intrusive_list_column& c2) { + if (&c1 == &c2) return true; + + if constexpr (Master_matrix::Option_list::is_z2) { + return c1.column_ == c2.column_; + } else { + auto it1 = c1.column_.begin(); + auto it2 = c2.column_.begin(); + if (c1.column_.size() != c2.column_.size()) return false; + while (it1 != c1.column_.end() && it2 != c2.column_.end()) { + if (it1->get_row_index() != it2->get_row_index() || it1->get_element() != it2->get_element()) return false; + ++it1; + ++it2; + } + return true; + } + } + friend bool operator<(const Intrusive_list_column& c1, const Intrusive_list_column& c2) { + if (&c1 == &c2) return false; + + if constexpr (Master_matrix::Option_list::is_z2) { + return c1.column_ < c2.column_; + } else { + auto it1 = c1.column_.begin(); + auto it2 = c2.column_.begin(); + while (it1 != c1.column_.end() && it2 != c2.column_.end()) { + if (it1->get_row_index() != it2->get_row_index()) return it1->get_row_index() < it2->get_row_index(); + if (it1->get_element() != it2->get_element()) return it1->get_element() < it2->get_element(); + ++it1; + ++it2; + } + return it2 != c2.column_.end(); + } + } + + // void set_operators(Field_operators* operators){ operators_ = operators; } + + // Disabled with row access. + Intrusive_list_column& operator=(const Intrusive_list_column& other); + + friend void swap(Intrusive_list_column& col1, Intrusive_list_column& col2) { + swap(static_cast(col1), + static_cast(col2)); + swap(static_cast(col1), + static_cast(col2)); + swap(static_cast(col1), + static_cast(col2)); + col1.column_.swap(col2.column_); + std::swap(col1.operators_, col2.operators_); + std::swap(col1.cellPool_, col2.cellPool_); + } + + private: + using ra_opt = typename Master_matrix::Row_access_option; + using dim_opt = typename Master_matrix::Column_dimension_option; + using chain_opt = typename Master_matrix::Chain_column_option; + + // Cloner object function for boost intrusive container + struct new_cloner { + new_cloner(Cell_constructor* cellPool) : cellPool_(cellPool){}; + + Cell* operator()(const Cell& clone_this) { return cellPool_->construct(clone_this); } + + Cell_constructor* cellPool_; + }; + + // The disposer object function for boost intrusive container + struct delete_disposer { + delete_disposer(){}; + delete_disposer(Intrusive_list_column* col) : col_(col){}; + + void operator()(Cell* delete_this) { + if constexpr (Master_matrix::Option_list::has_row_access) col_->unlink(delete_this); + col_->cellPool_->destroy(delete_this); + } + + Intrusive_list_column* col_; + }; + + Field_operators* operators_; + Cell_constructor* cellPool_; + Column_type column_; + + void _delete_cell(iterator& it); + void _insert_cell(const Field_element_type& value, id_index rowIndex, const iterator& position); + void _insert_cell(id_index rowIndex, const iterator& position); + template + bool _add(const Cell_range& column); + template + bool _multiply_and_add(const Field_element_type& val, const Cell_range& column); + template + bool _multiply_and_add(const Cell_range& column, const Field_element_type& val); + + void _verifyCellConstructor() { + if (cellPool_ == nullptr) { + if constexpr (std::is_same_v >) { + cellPool_ = &Master_matrix::defaultCellConstructor; + } else { + throw std::invalid_argument("Cell constructor pointer cannot be null."); + } + } + } }; -template -inline Intrusive_list_column::Intrusive_list_column(Field_operators* operators, Cell_constructor* cellConstructor) - : ra_opt(), dim_opt(), chain_opt(), operators_(operators), - cellPool_(cellConstructor), column_() -{ - if (operators_ == nullptr && cellPool_ == nullptr) return; //to allow default constructor which gives a dummy column - _verifyCellConstructor(); -} - -template -template -inline Intrusive_list_column::Intrusive_list_column(const Container_type &nonZeroRowIndices, Field_operators* operators, Cell_constructor* cellConstructor) - : ra_opt(), - dim_opt(nonZeroRowIndices.size() == 0 ? 0 : nonZeroRowIndices.size() - 1), - chain_opt(), - operators_(operators), - cellPool_(cellConstructor), - column_() -{ - static_assert(!Master_matrix::isNonBasic || Master_matrix::Option_list::is_of_boundary_type, - "Constructor not available for chain columns, please specify the dimension of the chain."); - - _verifyCellConstructor(); - - if constexpr (Master_matrix::Option_list::is_z2){ - for (id_index id : nonZeroRowIndices){ - _insert_cell(id, column_.end()); - } - } else { - for (const auto& p : nonZeroRowIndices){ - _insert_cell(operators_->get_value(p.second), p.first, column_.end()); - } - } -} - -template -template -inline Intrusive_list_column::Intrusive_list_column( - index columnIndex, const Container_type &nonZeroRowIndices, Row_container_type *rowContainer, Field_operators* operators, Cell_constructor* cellConstructor) - : ra_opt(columnIndex, rowContainer), - dim_opt(nonZeroRowIndices.size() == 0 ? 0 : nonZeroRowIndices.size() - 1), - chain_opt([&]{ - if constexpr (Master_matrix::Option_list::is_z2){ - return nonZeroRowIndices.begin() == nonZeroRowIndices.end() ? -1 : *std::prev(nonZeroRowIndices.end()); - } else { - return nonZeroRowIndices.begin() == nonZeroRowIndices.end() ? -1 : std::prev(nonZeroRowIndices.end())->first; - } - }()), - operators_(operators), - cellPool_(cellConstructor), - column_() -{ - static_assert(!Master_matrix::isNonBasic || Master_matrix::Option_list::is_of_boundary_type, - "Constructor not available for chain columns, please specify the dimension of the chain."); - - _verifyCellConstructor(); - - if constexpr (Master_matrix::Option_list::is_z2){ - for (id_index id : nonZeroRowIndices){ - _insert_cell(id, column_.end()); - } - } else { - for (const auto& p : nonZeroRowIndices){ - _insert_cell(operators_->get_value(p.second), p.first, column_.end()); - } - } -} - -template -template -inline Intrusive_list_column::Intrusive_list_column( - const Container_type &nonZeroRowIndices, dimension_type dimension, Field_operators* operators, Cell_constructor* cellConstructor) - : ra_opt(), - dim_opt(dimension), - chain_opt([&]{ - if constexpr (Master_matrix::Option_list::is_z2){ - return nonZeroRowIndices.begin() == nonZeroRowIndices.end() ? -1 : *std::prev(nonZeroRowIndices.end()); - } else { - return nonZeroRowIndices.begin() == nonZeroRowIndices.end() ? -1 : std::prev(nonZeroRowIndices.end())->first; - } - }()), - operators_(operators), - cellPool_(cellConstructor), - column_() -{ - _verifyCellConstructor(); - - if constexpr (Master_matrix::Option_list::is_z2){ - for (id_index id : nonZeroRowIndices){ - _insert_cell(id, column_.end()); - } - } else { - for (const auto& p : nonZeroRowIndices){ - _insert_cell(operators_->get_value(p.second), p.first, column_.end()); - } - } -} - -template -template -inline Intrusive_list_column::Intrusive_list_column( - index columnIndex, const Container_type &nonZeroRowIndices, dimension_type dimension, Row_container_type *rowContainer, Field_operators* operators, Cell_constructor* cellConstructor) - : ra_opt(columnIndex, rowContainer), - dim_opt(dimension), - chain_opt([&]{ - if constexpr (Master_matrix::Option_list::is_z2){ - return nonZeroRowIndices.begin() == nonZeroRowIndices.end() ? -1 : *std::prev(nonZeroRowIndices.end()); - } else { - return nonZeroRowIndices.begin() == nonZeroRowIndices.end() ? -1 : std::prev(nonZeroRowIndices.end())->first; - } - }()), - operators_(operators), - cellPool_(cellConstructor), - column_() -{ - _verifyCellConstructor(); - - if constexpr (Master_matrix::Option_list::is_z2){ - for (id_index id : nonZeroRowIndices){ - _insert_cell(id, column_.end()); - } - } else { - for (const auto& p : nonZeroRowIndices){ - _insert_cell(operators_->get_value(p.second), p.first, column_.end()); - } - } -} - -template -inline Intrusive_list_column::Intrusive_list_column(const Intrusive_list_column &column, Field_operators* operators, Cell_constructor* cellConstructor) - : ra_opt(), - dim_opt(static_cast(column)), - chain_opt(static_cast(column)), - operators_(operators == nullptr ? column.operators_ : operators), - cellPool_(cellConstructor == nullptr ? column.cellPool_ : cellConstructor), - column_() -{ - static_assert(!Master_matrix::Option_list::has_row_access, - "Simple copy constructor not available when row access option enabled. Please specify the new column index and the row container."); - column_.clone_from(column.column_, new_cloner(cellPool_), delete_disposer(this)); -} - -template -template -inline Intrusive_list_column::Intrusive_list_column( - const Intrusive_list_column &column, index columnIndex, Row_container_type *rowContainer, Field_operators* operators, Cell_constructor* cellConstructor) - : ra_opt(columnIndex, rowContainer), - dim_opt(static_cast(column)), - chain_opt(static_cast(column)), - operators_(operators == nullptr ? column.operators_ : operators), - cellPool_(cellConstructor == nullptr ? column.cellPool_ : cellConstructor), - column_() -{ - for (const Cell& cell : column.column_){ - if constexpr (Master_matrix::Option_list::is_z2){ - _insert_cell(cell.get_row_index(), column_.end()); - } else { - _insert_cell(cell.get_element(), cell.get_row_index(), column_.end()); - } - } -} - -template -inline Intrusive_list_column::Intrusive_list_column(Intrusive_list_column &&column) noexcept - : ra_opt(std::move(static_cast(column))), - dim_opt(std::move(static_cast(column))), - chain_opt(std::move(static_cast(column))), - operators_(std::exchange(column.operators_, nullptr)), - cellPool_(std::exchange(column.cellPool_, nullptr)), - column_(std::move(column.column_)) -{} +template +inline Intrusive_list_column::Intrusive_list_column(Field_operators* operators, + Cell_constructor* cellConstructor) + : ra_opt(), dim_opt(), chain_opt(), operators_(operators), cellPool_(cellConstructor), column_() +{ + if (operators_ == nullptr && cellPool_ == nullptr) return; //to allow default constructor which gives a dummy column + _verifyCellConstructor(); +} + +template +template +inline Intrusive_list_column::Intrusive_list_column( + const Container_type& nonZeroRowIndices, Field_operators* operators, Cell_constructor* cellConstructor) + : ra_opt(), + dim_opt(nonZeroRowIndices.size() == 0 ? 0 : nonZeroRowIndices.size() - 1), + chain_opt(), + operators_(operators), + cellPool_(cellConstructor), + column_() +{ + static_assert(!Master_matrix::isNonBasic || Master_matrix::Option_list::is_of_boundary_type, + "Constructor not available for chain columns, please specify the dimension of the chain."); + + _verifyCellConstructor(); + + if constexpr (Master_matrix::Option_list::is_z2) { + for (id_index id : nonZeroRowIndices) { + _insert_cell(id, column_.end()); + } + } else { + for (const auto& p : nonZeroRowIndices) { + _insert_cell(operators_->get_value(p.second), p.first, column_.end()); + } + } +} + +template +template +inline Intrusive_list_column::Intrusive_list_column( + index columnIndex, + const Container_type& nonZeroRowIndices, + Row_container_type* rowContainer, + Field_operators* operators, + Cell_constructor* cellConstructor) + : ra_opt(columnIndex, rowContainer), + dim_opt(nonZeroRowIndices.size() == 0 ? 0 : nonZeroRowIndices.size() - 1), + chain_opt([&] { + if constexpr (Master_matrix::Option_list::is_z2) { + return nonZeroRowIndices.begin() == nonZeroRowIndices.end() ? -1 : *std::prev(nonZeroRowIndices.end()); + } else { + return nonZeroRowIndices.begin() == nonZeroRowIndices.end() ? -1 : std::prev(nonZeroRowIndices.end())->first; + } + }()), + operators_(operators), + cellPool_(cellConstructor), + column_() +{ + static_assert(!Master_matrix::isNonBasic || Master_matrix::Option_list::is_of_boundary_type, + "Constructor not available for chain columns, please specify the dimension of the chain."); + + _verifyCellConstructor(); + + if constexpr (Master_matrix::Option_list::is_z2) { + for (id_index id : nonZeroRowIndices) { + _insert_cell(id, column_.end()); + } + } else { + for (const auto& p : nonZeroRowIndices) { + _insert_cell(operators_->get_value(p.second), p.first, column_.end()); + } + } +} + +template +template +inline Intrusive_list_column::Intrusive_list_column( + const Container_type& nonZeroRowIndices, + dimension_type dimension, + Field_operators* operators, + Cell_constructor* cellConstructor) + : ra_opt(), + dim_opt(dimension), + chain_opt([&] { + if constexpr (Master_matrix::Option_list::is_z2) { + return nonZeroRowIndices.begin() == nonZeroRowIndices.end() ? -1 : *std::prev(nonZeroRowIndices.end()); + } else { + return nonZeroRowIndices.begin() == nonZeroRowIndices.end() ? -1 : std::prev(nonZeroRowIndices.end())->first; + } + }()), + operators_(operators), + cellPool_(cellConstructor), + column_() +{ + _verifyCellConstructor(); + + if constexpr (Master_matrix::Option_list::is_z2) { + for (id_index id : nonZeroRowIndices) { + _insert_cell(id, column_.end()); + } + } else { + for (const auto& p : nonZeroRowIndices) { + _insert_cell(operators_->get_value(p.second), p.first, column_.end()); + } + } +} -template -inline Intrusive_list_column::~Intrusive_list_column() +template +template +inline Intrusive_list_column::Intrusive_list_column( + index columnIndex, + const Container_type& nonZeroRowIndices, + dimension_type dimension, + Row_container_type* rowContainer, + Field_operators* operators, + Cell_constructor* cellConstructor) + : ra_opt(columnIndex, rowContainer), + dim_opt(dimension), + chain_opt([&] { + if constexpr (Master_matrix::Option_list::is_z2) { + return nonZeroRowIndices.begin() == nonZeroRowIndices.end() ? -1 : *std::prev(nonZeroRowIndices.end()); + } else { + return nonZeroRowIndices.begin() == nonZeroRowIndices.end() ? -1 : std::prev(nonZeroRowIndices.end())->first; + } + }()), + operators_(operators), + cellPool_(cellConstructor), + column_() { - column_.clear_and_dispose(delete_disposer(this)); + _verifyCellConstructor(); + + if constexpr (Master_matrix::Option_list::is_z2) { + for (id_index id : nonZeroRowIndices) { + _insert_cell(id, column_.end()); + } + } else { + for (const auto& p : nonZeroRowIndices) { + _insert_cell(operators_->get_value(p.second), p.first, column_.end()); + } + } } -template -inline std::vector::Field_element_type> -Intrusive_list_column::get_content(int columnLength) const +template +inline Intrusive_list_column::Intrusive_list_column( + const Intrusive_list_column& column, Field_operators* operators, Cell_constructor* cellConstructor) + : ra_opt(), + dim_opt(static_cast(column)), + chain_opt(static_cast(column)), + operators_(operators == nullptr ? column.operators_ : operators), + cellPool_(cellConstructor == nullptr ? column.cellPool_ : cellConstructor), + column_() { - if (columnLength < 0 && column_.size() > 0) columnLength = column_.back().get_row_index() + 1; - else if (columnLength < 0) return std::vector(); + static_assert(!Master_matrix::Option_list::has_row_access, + "Simple copy constructor not available when row access option enabled. Please specify the new column " + "index and the row container."); + column_.clone_from(column.column_, new_cloner(cellPool_), delete_disposer(this)); +} - std::vector container(columnLength); - for (auto it = column_.begin(); it != column_.end() && it->get_row_index() < static_cast(columnLength); ++it){ - if constexpr (Master_matrix::Option_list::is_z2){ - container[it->get_row_index()] = 1; - } else { - container[it->get_row_index()] = it->get_element(); - } - } - return container; +template +template +inline Intrusive_list_column::Intrusive_list_column( + const Intrusive_list_column& column, index columnIndex, Row_container_type* rowContainer, + Field_operators* operators, Cell_constructor* cellConstructor) + : ra_opt(columnIndex, rowContainer), + dim_opt(static_cast(column)), + chain_opt(static_cast(column)), + operators_(operators == nullptr ? column.operators_ : operators), + cellPool_(cellConstructor == nullptr ? column.cellPool_ : cellConstructor), + column_() +{ + for (const Cell& cell : column.column_) { + if constexpr (Master_matrix::Option_list::is_z2) { + _insert_cell(cell.get_row_index(), column_.end()); + } else { + _insert_cell(cell.get_element(), cell.get_row_index(), column_.end()); + } + } } -template -inline bool Intrusive_list_column::is_non_zero(id_index rowIndex) const +template +inline Intrusive_list_column::Intrusive_list_column( + Intrusive_list_column&& column) noexcept + : ra_opt(std::move(static_cast(column))), + dim_opt(std::move(static_cast(column))), + chain_opt(std::move(static_cast(column))), + operators_(std::exchange(column.operators_, nullptr)), + cellPool_(std::exchange(column.cellPool_, nullptr)), + column_(std::move(column.column_)) +{} + +template +inline Intrusive_list_column::~Intrusive_list_column() { - //could be changed to dichotomic search as column is ordered by row index, - //but I am not sure if it is really worth it as there is no random access - //and the columns should not be that long anyway. - for (const Cell& cell : column_) - if (cell.get_row_index() == rowIndex) return true; + column_.clear_and_dispose(delete_disposer(this)); +} - return false; +template +inline std::vector::Field_element_type> +Intrusive_list_column::get_content(int columnLength) const +{ + if (columnLength < 0 && column_.size() > 0) + columnLength = column_.back().get_row_index() + 1; + else if (columnLength < 0) + return std::vector(); + + std::vector container(columnLength); + for (auto it = column_.begin(); it != column_.end() && it->get_row_index() < static_cast(columnLength); + ++it) { + if constexpr (Master_matrix::Option_list::is_z2) { + container[it->get_row_index()] = 1; + } else { + container[it->get_row_index()] = it->get_element(); + } + } + return container; +} + +template +inline bool Intrusive_list_column::is_non_zero(id_index rowIndex) const +{ + // could be changed to dichotomic search as column is ordered by row index, + // but I am not sure if it is really worth it as there is no random access + // and the columns should not be that long anyway. + for (const Cell& cell : column_) + if (cell.get_row_index() == rowIndex) return true; + + return false; +} + +template +inline bool Intrusive_list_column::is_empty() const +{ + return column_.empty(); } -template -inline bool Intrusive_list_column::is_empty() const +template +inline std::size_t Intrusive_list_column::size() const { - return column_.empty(); + return column_.size(); } -template -inline std::size_t Intrusive_list_column::size() const{ - return column_.size(); +template +template +inline void Intrusive_list_column::reorder(const Map_type& valueMap, + [[maybe_unused]] index columnIndex) +{ + static_assert(!Master_matrix::isNonBasic || Master_matrix::Option_list::is_of_boundary_type, + "Method not available for chain columns."); + + for (auto it = column_.begin(); it != column_.end(); ++it) { + Cell* cell = &(*it); + if constexpr (Master_matrix::Option_list::has_row_access) { + ra_opt::unlink(cell); + if (columnIndex != static_cast(-1)) cell->set_column_index(columnIndex); + } + cell->set_row_index(valueMap.at(cell->get_row_index())); + if constexpr (Master_matrix::Option_list::has_intrusive_rows && Master_matrix::Option_list::has_row_access) + ra_opt::insert_cell(cell->get_row_index(), cell); + } + + // all cells have to be deleted first, to avoid problem with insertion when row is a set + if constexpr (!Master_matrix::Option_list::has_intrusive_rows && Master_matrix::Option_list::has_row_access) { + for (auto it = column_.begin(); it != column_.end(); ++it) { + Cell* cell = &(*it); + ra_opt::insert_cell(cell->get_row_index(), cell); + } + } + + column_.sort(); } -template -template -inline void Intrusive_list_column::reorder(const Map_type &valueMap, [[maybe_unused]] index columnIndex) +template +inline void Intrusive_list_column::clear() { - static_assert(!Master_matrix::isNonBasic || Master_matrix::Option_list::is_of_boundary_type, - "Method not available for chain columns."); + static_assert(!Master_matrix::isNonBasic || Master_matrix::Option_list::is_of_boundary_type, + "Method not available for chain columns as a base element should not be empty."); - for (auto it = column_.begin(); it != column_.end(); ++it) { - Cell* cell = &(*it); - if constexpr (Master_matrix::Option_list::has_row_access) { - ra_opt::unlink(cell); - if (columnIndex != static_cast(-1)) cell->set_column_index(columnIndex); - } - cell->set_row_index(valueMap.at(cell->get_row_index())); - if constexpr (Master_matrix::Option_list::has_intrusive_rows && Master_matrix::Option_list::has_row_access) - ra_opt::insert_cell(cell->get_row_index(), cell); - } + column_.clear_and_dispose(delete_disposer(this)); +} - //all cells have to be deleted first, to avoid problem with insertion when row is a set - if constexpr (!Master_matrix::Option_list::has_intrusive_rows && Master_matrix::Option_list::has_row_access){ - for (auto it = column_.begin(); it != column_.end(); ++it) { - Cell* cell = &(*it); - ra_opt::insert_cell(cell->get_row_index(), cell); - } - } +template +inline void Intrusive_list_column::clear(id_index rowIndex) +{ + static_assert(!Master_matrix::isNonBasic || Master_matrix::Option_list::is_of_boundary_type, + "Method not available for chain columns."); - column_.sort(); + auto it = column_.begin(); + while (it != column_.end() && it->get_row_index() != rowIndex) it++; + if (it != column_.end()) _delete_cell(it); } -template -inline void Intrusive_list_column::clear() +template +inline typename Intrusive_list_column::id_index +Intrusive_list_column::get_pivot() const { - static_assert(!Master_matrix::isNonBasic || Master_matrix::Option_list::is_of_boundary_type, - "Method not available for chain columns as a base element should not be empty."); + static_assert(Master_matrix::isNonBasic, + "Method not available for base columns."); // could technically be, but is the notion usefull then? + + if constexpr (Master_matrix::Option_list::is_of_boundary_type) { + if (column_.empty()) return -1; + return column_.back().get_row_index(); + } else { + return chain_opt::get_pivot(); + } +} - column_.clear_and_dispose(delete_disposer(this)); +template +inline typename Intrusive_list_column::Field_element_type +Intrusive_list_column::get_pivot_value() const +{ + static_assert(Master_matrix::isNonBasic, + "Method not available for base columns."); // could technically be, but is the notion usefull then? + + if constexpr (Master_matrix::Option_list::is_z2) { + return 1; + } else { + if constexpr (Master_matrix::Option_list::is_of_boundary_type) { + if (column_.empty()) return 0; + return column_.back().get_element(); + } else { + if (chain_opt::get_pivot() == -1) return Field_element_type(); + for (const Cell& cell : column_) { + if (cell.get_row_index() == chain_opt::get_pivot()) return cell.get_element(); + } + return Field_element_type(); // should never happen if chain column is used properly + } + } } -template -inline void Intrusive_list_column::clear(id_index rowIndex) +template +inline typename Intrusive_list_column::iterator +Intrusive_list_column::begin() noexcept { - static_assert(!Master_matrix::isNonBasic || Master_matrix::Option_list::is_of_boundary_type, - "Method not available for chain columns."); + return column_.begin(); +} - auto it = column_.begin(); - while (it != column_.end() && it->get_row_index() != rowIndex) it++; - if (it != column_.end()) _delete_cell(it); +template +inline typename Intrusive_list_column::const_iterator +Intrusive_list_column::begin() const noexcept +{ + return column_.begin(); } -template -inline typename Intrusive_list_column::id_index Intrusive_list_column::get_pivot() const +template +inline typename Intrusive_list_column::iterator +Intrusive_list_column::end() noexcept { - static_assert(Master_matrix::isNonBasic, "Method not available for base columns."); //could technically be, but is the notion usefull then? + return column_.end(); +} - if constexpr (Master_matrix::Option_list::is_of_boundary_type){ - if (column_.empty()) return -1; - return column_.back().get_row_index(); - } else { - return chain_opt::get_pivot(); - } +template +inline typename Intrusive_list_column::const_iterator +Intrusive_list_column::end() const noexcept +{ + return column_.end(); } -template -inline typename Intrusive_list_column::Field_element_type Intrusive_list_column::get_pivot_value() const +template +inline typename Intrusive_list_column::reverse_iterator +Intrusive_list_column::rbegin() noexcept { - static_assert(Master_matrix::isNonBasic, "Method not available for base columns."); //could technically be, but is the notion usefull then? + return column_.rbegin(); +} - if constexpr (Master_matrix::Option_list::is_z2){ - return 1; - } else { - if constexpr (Master_matrix::Option_list::is_of_boundary_type){ - if (column_.empty()) return 0; - return column_.back().get_element(); - } else { - if (chain_opt::get_pivot() == -1) return Field_element_type(); - for (const Cell& cell : column_){ - if (cell.get_row_index() == chain_opt::get_pivot()) return cell.get_element(); - } - return Field_element_type(); //should never happen if chain column is used properly - } - } +template +inline typename Intrusive_list_column::const_reverse_iterator +Intrusive_list_column::rbegin() const noexcept +{ + return column_.rbegin(); } -template -inline typename Intrusive_list_column::iterator -Intrusive_list_column::begin() noexcept +template +inline typename Intrusive_list_column::reverse_iterator +Intrusive_list_column::rend() noexcept { - return column_.begin(); + return column_.rend(); } -template -inline typename Intrusive_list_column::const_iterator -Intrusive_list_column::begin() const noexcept +template +inline typename Intrusive_list_column::const_reverse_iterator +Intrusive_list_column::rend() const noexcept { - return column_.begin(); + return column_.rend(); +} + +template +template +inline Intrusive_list_column& +Intrusive_list_column::operator+=(const Cell_range& column) +{ + static_assert((!Master_matrix::isNonBasic || std::is_same_v), + "For boundary columns, the range has to be a column of same type to help ensure the validity of the " + "base element."); // could be removed, if we give the responsability to the user. + static_assert((!Master_matrix::isNonBasic || Master_matrix::Option_list::is_of_boundary_type), + "For chain columns, the given column cannot be constant."); + + _add(column); + + return *this; } -template -inline typename Intrusive_list_column::iterator -Intrusive_list_column::end() noexcept +template +inline Intrusive_list_column& +Intrusive_list_column::operator+=(Intrusive_list_column& column) { - return column_.end(); + if constexpr (Master_matrix::isNonBasic && !Master_matrix::Option_list::is_of_boundary_type) { + // assumes that the addition never zeros out this column. + if (_add(column)) { + chain_opt::swap_pivots(column); + dim_opt::swap_dimension(column); + } + } else { + _add(column); + } + + return *this; } -template -inline typename Intrusive_list_column::const_iterator -Intrusive_list_column::end() const noexcept +template +inline Intrusive_list_column& +Intrusive_list_column::operator*=(const Field_element_type& val) { - return column_.end(); + if constexpr (Master_matrix::Option_list::is_z2) { + if (val % 2 == 0) { + if constexpr (Master_matrix::isNonBasic && !Master_matrix::Option_list::is_of_boundary_type) { + throw std::invalid_argument("A chain column should not be multiplied by 0."); + } else { + clear(); + } + } + } else { + Field_element_type realVal = operators_->get_value(val); + + if (realVal == Field_operators::get_additive_identity()) { + if constexpr (Master_matrix::isNonBasic && !Master_matrix::Option_list::is_of_boundary_type) { + throw std::invalid_argument("A chain column should not be multiplied by 0."); + } else { + clear(); + } + return *this; + } + + if (realVal == Field_operators::get_multiplicative_identity()) return *this; + + for (Cell& cell : column_) { + cell.get_element() = operators_->multiply(cell.get_element(), realVal); + if constexpr (Master_matrix::Option_list::has_row_access) ra_opt::update_cell(cell); + } + } + + return *this; } -template -inline typename Intrusive_list_column::reverse_iterator -Intrusive_list_column::rbegin() noexcept +template +template +inline Intrusive_list_column& +Intrusive_list_column::multiply_and_add(const Field_element_type& val, + const Cell_range& column) { - return column_.rbegin(); + static_assert((!Master_matrix::isNonBasic || std::is_same_v), + "For boundary columns, the range has to be a column of same type to help ensure the validity of the " + "base element."); // could be removed, if we give the responsability to the user. + static_assert((!Master_matrix::isNonBasic || Master_matrix::Option_list::is_of_boundary_type), + "For chain columns, the given column cannot be constant."); + + if constexpr (Master_matrix::Option_list::is_z2) { + if (val) { + _add(column); + } else { + clear(); + _add(column); + } + } else { + _multiply_and_add(val, column); + } + + return *this; } -template -inline typename Intrusive_list_column::const_reverse_iterator -Intrusive_list_column::rbegin() const noexcept +template +inline Intrusive_list_column& +Intrusive_list_column::multiply_and_add(const Field_element_type& val, + Intrusive_list_column& column) { - return column_.rbegin(); + if constexpr (Master_matrix::isNonBasic && !Master_matrix::Option_list::is_of_boundary_type) { + // assumes that the addition never zeros out this column. + if constexpr (Master_matrix::Option_list::is_z2) { + if (val) { + if (_add(column)) { + chain_opt::swap_pivots(column); + dim_opt::swap_dimension(column); + } + } else { + throw std::invalid_argument("A chain column should not be multiplied by 0."); + } + } else { + if (_multiply_and_add(val, column)) { + chain_opt::swap_pivots(column); + dim_opt::swap_dimension(column); + } + } + } else { + if constexpr (Master_matrix::Option_list::is_z2) { + if (val) { + _add(column); + } else { + clear(); + _add(column); + } + } else { + _multiply_and_add(val, column); + } + } + + return *this; } -template -inline typename Intrusive_list_column::reverse_iterator -Intrusive_list_column::rend() noexcept +template +template +inline Intrusive_list_column& +Intrusive_list_column::multiply_and_add(const Cell_range& column, + const Field_element_type& val) { - return column_.rend(); + static_assert((!Master_matrix::isNonBasic || std::is_same_v), + "For boundary columns, the range has to be a column of same type to help ensure the validity of the " + "base element."); // could be removed, if we give the responsability to the user. + static_assert((!Master_matrix::isNonBasic || Master_matrix::Option_list::is_of_boundary_type), + "For chain columns, the given column cannot be constant."); + + if constexpr (Master_matrix::Option_list::is_z2) { + if (val) { + _add(column); + } + } else { + _multiply_and_add(column, val); + } + + return *this; } -template -inline typename Intrusive_list_column::const_reverse_iterator -Intrusive_list_column::rend() const noexcept +template +inline Intrusive_list_column& +Intrusive_list_column::multiply_and_add(Intrusive_list_column& column, + const Field_element_type& val) { - return column_.rend(); + if constexpr (Master_matrix::isNonBasic && !Master_matrix::Option_list::is_of_boundary_type) { + // assumes that the addition never zeros out this column. + if constexpr (Master_matrix::Option_list::is_z2) { + if (val) { + if (_add(column)) { + chain_opt::swap_pivots(column); + dim_opt::swap_dimension(column); + } + } + } else { + if (_multiply_and_add(column, val)) { + chain_opt::swap_pivots(column); + dim_opt::swap_dimension(column); + } + } + } else { + if constexpr (Master_matrix::Option_list::is_z2) { + if (val) { + _add(column); + } + } else { + _multiply_and_add(column, val); + } + } + + return *this; } -template -template -inline Intrusive_list_column & -Intrusive_list_column::operator+=(const Cell_range &column) +template +inline Intrusive_list_column& +Intrusive_list_column::operator=(const Intrusive_list_column& other) { - static_assert((!Master_matrix::isNonBasic || std::is_same_v), - "For boundary columns, the range has to be a column of same type to help ensure the validity of the base element."); //could be removed, if we give the responsability to the user. - static_assert((!Master_matrix::isNonBasic || Master_matrix::Option_list::is_of_boundary_type), - "For chain columns, the given column cannot be constant."); + static_assert(!Master_matrix::Option_list::has_row_access, "= assignement not enabled with row access option."); + + dim_opt::operator=(other); + chain_opt::operator=(other); - _add(column); + // order is important + column_.clear_and_dispose(delete_disposer(this)); + operators_ = other.operators_; + cellPool_ = other.cellPool_; + column_.clone_from(other.column_, new_cloner(cellPool_), delete_disposer(this)); - return *this; + return *this; } -template -inline Intrusive_list_column & -Intrusive_list_column::operator+=(Intrusive_list_column &column) +template +inline void Intrusive_list_column::_delete_cell(iterator& it) { - if constexpr (Master_matrix::isNonBasic && !Master_matrix::Option_list::is_of_boundary_type){ - //assumes that the addition never zeros out this column. - if (_add(column)) { - chain_opt::swap_pivots(column); - dim_opt::swap_dimension(column); - } - } else { - _add(column); - } + it = column_.erase_and_dispose(it, delete_disposer(this)); +} - return *this; +template +inline void Intrusive_list_column::_insert_cell(const Field_element_type& value, + id_index rowIndex, + const iterator& position) +{ + if constexpr (Master_matrix::Option_list::has_row_access) { + Cell* new_cell = cellPool_->construct(ra_opt::columnIndex_, rowIndex); + new_cell->set_element(value); + column_.insert(position, *new_cell); + ra_opt::insert_cell(rowIndex, new_cell); + } else { + Cell* new_cell = cellPool_->construct(rowIndex); + new_cell->set_element(value); + column_.insert(position, *new_cell); + } } -template -inline Intrusive_list_column& -Intrusive_list_column::operator*=(const Field_element_type& val) +template +inline void Intrusive_list_column::_insert_cell(id_index rowIndex, + const iterator& position) { - if constexpr (Master_matrix::Option_list::is_z2){ - if (val % 2 == 0){ - if constexpr (Master_matrix::isNonBasic && !Master_matrix::Option_list::is_of_boundary_type){ - throw std::invalid_argument("A chain column should not be multiplied by 0."); - } else { - clear(); - } - } - } else { - Field_element_type realVal = operators_->get_value(val); + if constexpr (Master_matrix::Option_list::has_row_access) { + Cell* new_cell = cellPool_->construct(ra_opt::columnIndex_, rowIndex); + column_.insert(position, *new_cell); + ra_opt::insert_cell(rowIndex, new_cell); + } else { + Cell* new_cell = cellPool_->construct(rowIndex); + column_.insert(position, *new_cell); + } +} - if (realVal == Field_operators::get_additive_identity()) { - if constexpr (Master_matrix::isNonBasic && !Master_matrix::Option_list::is_of_boundary_type){ - throw std::invalid_argument("A chain column should not be multiplied by 0."); - } else { - clear(); - } - return *this; - } +template +template +inline bool Intrusive_list_column::_add(const Cell_range& column) +{ + auto itTarget = column_.begin(); + auto itSource = column.begin(); + bool pivotIsZeroed = false; + + while (itTarget != column_.end() && itSource != column.end()) { + if (itTarget->get_row_index() < itSource->get_row_index()) { + ++itTarget; + } else if (itTarget->get_row_index() > itSource->get_row_index()) { + if constexpr (Master_matrix::Option_list::is_z2) { + _insert_cell(itSource->get_row_index(), itTarget); + } else { + _insert_cell(itSource->get_element(), itSource->get_row_index(), itTarget); + } + ++itSource; + } else { + if constexpr (Master_matrix::Option_list::is_z2) { + if constexpr (Master_matrix::isNonBasic && !Master_matrix::Option_list::is_of_boundary_type) { + if (itTarget->get_row_index() == chain_opt::get_pivot()) pivotIsZeroed = true; + } + _delete_cell(itTarget); + } else { + itTarget->get_element() = operators_->add(itTarget->get_element(), itSource->get_element()); + if (itTarget->get_element() == + Field_operators::get_additive_identity()) { // get_element is already modulo, so '==' works. + if constexpr (Master_matrix::isNonBasic && !Master_matrix::Option_list::is_of_boundary_type) { + if (itTarget->get_row_index() == chain_opt::get_pivot()) pivotIsZeroed = true; + } + _delete_cell(itTarget); + } else { + if constexpr (Master_matrix::Option_list::has_row_access) ra_opt::update_cell(*itTarget); + ++itTarget; + } + } + ++itSource; + } + } + + while (itSource != column.end()) { + if constexpr (Master_matrix::Option_list::is_z2) { + _insert_cell(itSource->get_row_index(), column_.end()); + } else { + _insert_cell(itSource->get_element(), itSource->get_row_index(), column_.end()); + } + ++itSource; + } + + return pivotIsZeroed; +} - if (realVal == Field_operators::get_multiplicative_identity()) return *this; +template +template +inline bool Intrusive_list_column::_multiply_and_add(const Field_element_type& val, + const Cell_range& column) +{ + bool pivotIsZeroed = false; + + if (val == 0u) { + if constexpr (Master_matrix::isNonBasic && !Master_matrix::Option_list::is_of_boundary_type) { + throw std::invalid_argument("A chain column should not be multiplied by 0."); + // this would not only mess up the base, but also the pivots stored. + } else { + clear(); + } + } + + auto itTarget = column_.begin(); + auto itSource = column.begin(); + while (itTarget != column_.end() && itSource != column.end()) { + if (itTarget->get_row_index() < itSource->get_row_index()) { + itTarget->get_element() = operators_->multiply(itTarget->get_element(), val); + if constexpr (Master_matrix::Option_list::has_row_access) ra_opt::update_cell(*itTarget); + ++itTarget; + } else if (itTarget->get_row_index() > itSource->get_row_index()) { + _insert_cell(itSource->get_element(), itSource->get_row_index(), itTarget); + ++itSource; + } else { + itTarget->get_element() = operators_->multiply_and_add(itTarget->get_element(), val, itSource->get_element()); + if (itTarget->get_element() == Field_operators::get_additive_identity()) { + if constexpr (Master_matrix::isNonBasic && !Master_matrix::Option_list::is_of_boundary_type) { + if (itTarget->get_row_index() == chain_opt::get_pivot()) pivotIsZeroed = true; + } + _delete_cell(itTarget); + } else { + if constexpr (Master_matrix::Option_list::has_row_access) ra_opt::update_cell(*itTarget); + ++itTarget; + } + ++itSource; + } + } + + while (itTarget != column_.end()) { + itTarget->get_element() = operators_->multiply(itTarget->get_element(), val); + if constexpr (Master_matrix::Option_list::has_row_access) ra_opt::update_cell(*itTarget); + itTarget++; + } + + while (itSource != column.end()) { + _insert_cell(itSource->get_element(), itSource->get_row_index(), column_.end()); + ++itSource; + } + + return pivotIsZeroed; +} + +template +template +inline bool Intrusive_list_column::_multiply_and_add(const Cell_range& column, + const Field_element_type& val) +{ + if (val == 0u) { + return false; + } + + bool pivotIsZeroed = false; + + auto itTarget = column_.begin(); + auto itSource = column.begin(); + while (itTarget != column_.end() && itSource != column.end()) { + if (itTarget->get_row_index() < itSource->get_row_index()) { + ++itTarget; + } else if (itTarget->get_row_index() > itSource->get_row_index()) { + _insert_cell(operators_->multiply(itSource->get_element(), val), itSource->get_row_index(), itTarget); + ++itSource; + } else { + itTarget->get_element() = operators_->multiply_and_add(itSource->get_element(), val, itTarget->get_element()); + if (itTarget->get_element() == Field_operators::get_additive_identity()) { + if constexpr (Master_matrix::isNonBasic && !Master_matrix::Option_list::is_of_boundary_type) { + if (itTarget->get_row_index() == chain_opt::get_pivot()) pivotIsZeroed = true; + } + _delete_cell(itTarget); + } else { + if constexpr (Master_matrix::Option_list::has_row_access) ra_opt::update_cell(*itTarget); + ++itTarget; + } + ++itSource; + } + } + + while (itSource != column.end()) { + _insert_cell(operators_->multiply(itSource->get_element(), val), itSource->get_row_index(), column_.end()); + ++itSource; + } + + return pivotIsZeroed; +} - for (Cell& cell : column_){ - cell.get_element() = operators_->multiply(cell.get_element(), realVal); - if constexpr (Master_matrix::Option_list::has_row_access) - ra_opt::update_cell(cell); - } - } - - return *this; -} - -template -template -inline Intrusive_list_column& -Intrusive_list_column::multiply_and_add(const Field_element_type& val, const Cell_range& column) -{ - static_assert((!Master_matrix::isNonBasic || std::is_same_v), - "For boundary columns, the range has to be a column of same type to help ensure the validity of the base element."); //could be removed, if we give the responsability to the user. - static_assert((!Master_matrix::isNonBasic || Master_matrix::Option_list::is_of_boundary_type), - "For chain columns, the given column cannot be constant."); - - if constexpr (Master_matrix::Option_list::is_z2){ - if (val){ - _add(column); - } else { - clear(); - _add(column); - } - } else { - _multiply_and_add(val, column); - } - - return *this; -} - -template -inline Intrusive_list_column & -Intrusive_list_column::multiply_and_add(const Field_element_type& val, Intrusive_list_column& column) -{ - if constexpr (Master_matrix::isNonBasic && !Master_matrix::Option_list::is_of_boundary_type){ - //assumes that the addition never zeros out this column. - if constexpr (Master_matrix::Option_list::is_z2){ - if (val){ - if (_add(column)){ - chain_opt::swap_pivots(column); - dim_opt::swap_dimension(column); - } - } else { - throw std::invalid_argument("A chain column should not be multiplied by 0."); - } - } else { - if (_multiply_and_add(val, column)){ - chain_opt::swap_pivots(column); - dim_opt::swap_dimension(column); - } - } - } else { - if constexpr (Master_matrix::Option_list::is_z2){ - if (val){ - _add(column); - } else { - clear(); - _add(column); - } - } else { - _multiply_and_add(val, column); - } - } - - return *this; -} - -template -template -inline Intrusive_list_column & -Intrusive_list_column::multiply_and_add(const Cell_range& column, const Field_element_type& val) -{ - static_assert((!Master_matrix::isNonBasic || std::is_same_v), - "For boundary columns, the range has to be a column of same type to help ensure the validity of the base element."); //could be removed, if we give the responsability to the user. - static_assert((!Master_matrix::isNonBasic || Master_matrix::Option_list::is_of_boundary_type), - "For chain columns, the given column cannot be constant."); - - if constexpr (Master_matrix::Option_list::is_z2){ - if (val){ - _add(column); - } - } else { - _multiply_and_add(column, val); - } - - return *this; -} - -template -inline Intrusive_list_column & -Intrusive_list_column::multiply_and_add(Intrusive_list_column& column, const Field_element_type& val) -{ - if constexpr (Master_matrix::isNonBasic && !Master_matrix::Option_list::is_of_boundary_type){ - //assumes that the addition never zeros out this column. - if constexpr (Master_matrix::Option_list::is_z2){ - if (val){ - if (_add(column)){ - chain_opt::swap_pivots(column); - dim_opt::swap_dimension(column); - } - } - } else { - if (_multiply_and_add(column, val)){ - chain_opt::swap_pivots(column); - dim_opt::swap_dimension(column); - } - } - } else { - if constexpr (Master_matrix::Option_list::is_z2){ - if (val){ - _add(column); - } - } else { - _multiply_and_add(column, val); - } - } - - return *this; -} - -template -inline Intrusive_list_column & -Intrusive_list_column::operator=(const Intrusive_list_column& other) -{ - static_assert(!Master_matrix::Option_list::has_row_access, "= assignement not enabled with row access option."); - - dim_opt::operator=(other); - chain_opt::operator=(other); - - //order is important - column_.clear_and_dispose(delete_disposer(this)); - operators_ = other.operators_; - cellPool_ = other.cellPool_; - column_.clone_from(other.column_, new_cloner(cellPool_), delete_disposer(this)); - - return *this; -} - -template -inline void Intrusive_list_column::_delete_cell(iterator &it) -{ - it = column_.erase_and_dispose(it, delete_disposer(this)); -} - -template -inline void Intrusive_list_column::_insert_cell( - const Field_element_type &value, id_index rowIndex, const iterator &position) -{ - if constexpr (Master_matrix::Option_list::has_row_access){ - Cell *new_cell = cellPool_->construct(ra_opt::columnIndex_, rowIndex); - new_cell->set_element(value); - column_.insert(position, *new_cell); - ra_opt::insert_cell(rowIndex, new_cell); - } else { - Cell *new_cell = cellPool_->construct(rowIndex); - new_cell->set_element(value); - column_.insert(position, *new_cell); - } -} - -template -inline void Intrusive_list_column::_insert_cell( - id_index rowIndex, const iterator &position) -{ - if constexpr (Master_matrix::Option_list::has_row_access){ - Cell *new_cell = cellPool_->construct(ra_opt::columnIndex_, rowIndex); - column_.insert(position, *new_cell); - ra_opt::insert_cell(rowIndex, new_cell); - } else { - Cell *new_cell = cellPool_->construct(rowIndex); - column_.insert(position, *new_cell); - } -} - -template -template -inline bool Intrusive_list_column::_add(const Cell_range &column) -{ - auto itTarget = column_.begin(); - auto itSource = column.begin(); - bool pivotIsZeroed = false; - - while (itTarget != column_.end() && itSource != column.end()) - { - if (itTarget->get_row_index() < itSource->get_row_index()) { - ++itTarget; - } else if (itTarget->get_row_index() > itSource->get_row_index()) { - if constexpr (Master_matrix::Option_list::is_z2){ - _insert_cell(itSource->get_row_index(), itTarget); - } else { - _insert_cell(itSource->get_element(), itSource->get_row_index(), itTarget); - } - ++itSource; - } else { - if constexpr (Master_matrix::Option_list::is_z2){ - if constexpr (Master_matrix::isNonBasic && !Master_matrix::Option_list::is_of_boundary_type){ - if (itTarget->get_row_index() == chain_opt::get_pivot()) pivotIsZeroed = true; - } - _delete_cell(itTarget); - } else { - itTarget->get_element() = operators_->add(itTarget->get_element(), itSource->get_element()); - if (itTarget->get_element() == Field_operators::get_additive_identity()){ //get_element is already modulo, so '==' works. - if constexpr (Master_matrix::isNonBasic && !Master_matrix::Option_list::is_of_boundary_type){ - if (itTarget->get_row_index() == chain_opt::get_pivot()) pivotIsZeroed = true; - } - _delete_cell(itTarget); - } else { - if constexpr (Master_matrix::Option_list::has_row_access) - ra_opt::update_cell(*itTarget); - ++itTarget; - } - } - ++itSource; - } - } - - while (itSource != column.end()) { - if constexpr (Master_matrix::Option_list::is_z2){ - _insert_cell(itSource->get_row_index(), column_.end()); - } else { - _insert_cell(itSource->get_element(), itSource->get_row_index(), column_.end()); - } - ++itSource; - } - - return pivotIsZeroed; -} - -template -template -inline bool Intrusive_list_column::_multiply_and_add(const Field_element_type& val, const Cell_range& column) -{ - bool pivotIsZeroed = false; - - if (val == 0u) { - if constexpr (Master_matrix::isNonBasic && !Master_matrix::Option_list::is_of_boundary_type){ - throw std::invalid_argument("A chain column should not be multiplied by 0."); - //this would not only mess up the base, but also the pivots stored. - } else { - clear(); - } - } - - auto itTarget = column_.begin(); - auto itSource = column.begin(); - while (itTarget != column_.end() && itSource != column.end()) - { - if (itTarget->get_row_index() < itSource->get_row_index()) { - itTarget->get_element() = operators_->multiply(itTarget->get_element(), val); - if constexpr (Master_matrix::Option_list::has_row_access) - ra_opt::update_cell(*itTarget); - ++itTarget; - } else if (itTarget->get_row_index() > itSource->get_row_index()) { - _insert_cell(itSource->get_element(), itSource->get_row_index(), itTarget); - ++itSource; - } else { - itTarget->get_element() = operators_->multiply_and_add(itTarget->get_element(), val, itSource->get_element()); - if (itTarget->get_element() == Field_operators::get_additive_identity()){ - if constexpr (Master_matrix::isNonBasic && !Master_matrix::Option_list::is_of_boundary_type){ - if (itTarget->get_row_index() == chain_opt::get_pivot()) pivotIsZeroed = true; - } - _delete_cell(itTarget); - } else { - if constexpr (Master_matrix::Option_list::has_row_access) - ra_opt::update_cell(*itTarget); - ++itTarget; - } - ++itSource; - } - } - - while (itTarget != column_.end()){ - itTarget->get_element() = operators_->multiply(itTarget->get_element(), val); - if constexpr (Master_matrix::Option_list::has_row_access) - ra_opt::update_cell(*itTarget); - itTarget++; - } - - while (itSource != column.end()) { - _insert_cell(itSource->get_element(), itSource->get_row_index(), column_.end()); - ++itSource; - } - - return pivotIsZeroed; -} - -template -template -inline bool Intrusive_list_column::_multiply_and_add(const Cell_range& column, const Field_element_type& val) -{ - if (val == 0u) { - return false; - } - - bool pivotIsZeroed = false; - - auto itTarget = column_.begin(); - auto itSource = column.begin(); - while (itTarget != column_.end() && itSource != column.end()) - { - if (itTarget->get_row_index() < itSource->get_row_index()) { - ++itTarget; - } else if (itTarget->get_row_index() > itSource->get_row_index()) { - _insert_cell(operators_->multiply(itSource->get_element(), val), itSource->get_row_index(), itTarget); - ++itSource; - } else { - itTarget->get_element() = operators_->multiply_and_add(itSource->get_element(), val, itTarget->get_element()); - if (itTarget->get_element() == Field_operators::get_additive_identity()){ - if constexpr (Master_matrix::isNonBasic && !Master_matrix::Option_list::is_of_boundary_type){ - if (itTarget->get_row_index() == chain_opt::get_pivot()) pivotIsZeroed = true; - } - _delete_cell(itTarget); - } else { - if constexpr (Master_matrix::Option_list::has_row_access) - ra_opt::update_cell(*itTarget); - ++itTarget; - } - ++itSource; - } - } - - while (itSource != column.end()) { - _insert_cell(operators_->multiply(itSource->get_element(), val), itSource->get_row_index(), column_.end()); - ++itSource; - } - - return pivotIsZeroed; -} - -} //namespace persistence_matrix -} //namespace Gudhi - -template -struct std::hash > -{ - size_t operator()(const Gudhi::persistence_matrix::Intrusive_list_column& column) const - { - std::size_t seed = 0; - for (auto& cell : column){ - seed ^= std::hash()(cell.get_row_index() * static_cast(cell.get_element())) + 0x9e3779b9 + (seed << 6) + (seed >> 2); - } - return seed; - } +} // namespace persistence_matrix +} // namespace Gudhi + +/** + * @brief Hash method for @ref Gudhi::persistence_matrix::Intrusive_list_column. + * + * @tparam Master_matrix Template parameter of @ref Gudhi::persistence_matrix::Intrusive_list_column. + * @tparam Cell_constructor Template parameter of @ref Gudhi::persistence_matrix::Intrusive_list_column. + */ +template +struct std::hash > +{ + size_t operator()( + const Gudhi::persistence_matrix::Intrusive_list_column& column) const { + std::size_t seed = 0; + for (auto& cell : column) { + seed ^= std::hash()(cell.get_row_index() * static_cast(cell.get_element())) + + 0x9e3779b9 + (seed << 6) + (seed >> 2); + } + return seed; + } }; -#endif // PM_INTRUSIVE_LIST_COLUMN_H +#endif // PM_INTRUSIVE_LIST_COLUMN_H diff --git a/src/Persistence_matrix/include/gudhi/Persistence_matrix/columns/intrusive_set_column.h b/src/Persistence_matrix/include/gudhi/Persistence_matrix/columns/intrusive_set_column.h index e3d05646b1..3fecff3456 100644 --- a/src/Persistence_matrix/include/gudhi/Persistence_matrix/columns/intrusive_set_column.h +++ b/src/Persistence_matrix/include/gudhi/Persistence_matrix/columns/intrusive_set_column.h @@ -2,19 +2,26 @@ * 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-24 Inria * * Modification(s): * - YYYY/MM Author: Description of the modification */ +/** + * @file intrusive_set_column.h + * @author Hannah Schreiber + * @brief Contains the @ref Intrusive_set_column class. + * Also defines the std::hash method for @ref Intrusive_set_column. + */ + #ifndef PM_INTRUSIVE_SET_COLUMN_H #define PM_INTRUSIVE_SET_COLUMN_H #include #include #include -#include //std::swap, std::move & std::exchange +#include //std::swap, std::move & std::exchange #include @@ -24,939 +31,992 @@ namespace Gudhi { namespace persistence_matrix { -template > -class Intrusive_set_column : public Master_matrix::Row_access_option, - public Master_matrix::Column_dimension_option, - public Master_matrix::Chain_column_option -{ -public: - using Master = Master_matrix; - using Field_operators = typename Master_matrix::Field_operators; - using Field_element_type = typename Master_matrix::element_type; - using index = typename Master_matrix::index; - using id_index = typename Master_matrix::id_index; - using dimension_type = typename Master_matrix::dimension_type; - - using Cell = typename Master_matrix::Cell_type; - using Column_type = boost::intrusive::set < - Cell, - boost::intrusive::constant_time_size, - boost::intrusive::base_hook< typename Master_matrix::base_hook_matrix_set_column > - >; - using iterator = typename Column_type::iterator; - using const_iterator = typename Column_type::const_iterator; - using reverse_iterator = typename Column_type::reverse_iterator; - using const_reverse_iterator = typename Column_type::const_reverse_iterator; - - Intrusive_set_column(Field_operators* operators = nullptr, Cell_constructor* cellConstructor = nullptr); - template - Intrusive_set_column(const Container_type& nonZeroRowIndices, Field_operators* operators, Cell_constructor* cellConstructor); //has to be a boundary for boundary, has no sense for chain if dimension is needed - template - Intrusive_set_column(index columnIndex, const Container_type& nonZeroRowIndices, Row_container_type *rowContainer, Field_operators* operators, Cell_constructor* cellConstructor); //has to be a boundary for boundary, has no sense for chain if dimension is needed - template - Intrusive_set_column(const Container_type& nonZeroChainRowIndices, dimension_type dimension, Field_operators* operators, Cell_constructor* cellConstructor); //dimension gets ignored for base - template - Intrusive_set_column(index columnIndex, const Container_type& nonZeroChainRowIndices, dimension_type dimension, Row_container_type *rowContainer, Field_operators* operators, Cell_constructor* cellConstructor); //dimension gets ignored for base - Intrusive_set_column(const Intrusive_set_column& column, Field_operators* operators = nullptr, Cell_constructor* cellConstructor = nullptr); - template - Intrusive_set_column(const Intrusive_set_column& column, index columnIndex, Row_container_type *rowContainer, Field_operators* operators = nullptr, Cell_constructor* cellConstructor = nullptr); - Intrusive_set_column(Intrusive_set_column&& column) noexcept; - ~Intrusive_set_column(); - - std::vector get_content(int columnLength = -1) const; - bool is_non_zero(id_index rowIndex) const; - bool is_empty() const; - std::size_t size() const; - - //**************** - //only for base and boundary - template - void reorder(const Map_type& valueMap, [[maybe_unused]] index columnIndex = -1); //used for lazy row swaps - void clear(); - void clear(id_index rowIndex); - //**************** - - //**************** - //only for chain and boundary - id_index get_pivot() const; - Field_element_type get_pivot_value() const; - //**************** - - iterator begin() noexcept; - const_iterator begin() const noexcept; - iterator end() noexcept; - const_iterator end() const noexcept; - reverse_iterator rbegin() noexcept; - const_reverse_iterator rbegin() const noexcept; - reverse_iterator rend() noexcept; - const_reverse_iterator rend() const noexcept; - - template - Intrusive_set_column& operator+=(const Cell_range& column); //for base & boundary except vector - Intrusive_set_column& operator+=(Intrusive_set_column &column); //for chain and vector - - Intrusive_set_column& operator*=(unsigned int v); - - //this = v * this + column - template - Intrusive_set_column& multiply_and_add(const Field_element_type& val, const Cell_range& column); //for base & boundary except vector - Intrusive_set_column& multiply_and_add(const Field_element_type& val, Intrusive_set_column& column); //for chain and vector - //this = this + column * v - template - Intrusive_set_column& multiply_and_add(const Cell_range& column, const Field_element_type& val); //for base & boundary except vector - Intrusive_set_column& multiply_and_add(Intrusive_set_column& column, const Field_element_type& val); //for chain and vector - - friend bool operator==(const Intrusive_set_column& c1, const Intrusive_set_column& c2){ - if (&c1 == &c2) return true; - - if constexpr (Master_matrix::Option_list::is_z2){ - return c1.column_ == c2.column_; - } else { - auto it1 = c1.column_.begin(); - auto it2 = c2.column_.begin(); - if (c1.column_.size() != c2.column_.size()) return false; - while (it1 != c1.column_.end() && it2 != c2.column_.end()) { - if (it1->get_row_index() != it2->get_row_index() || it1->get_element() != it2->get_element()) - return false; - ++it1; ++it2; - } - return true; - } - } - friend bool operator<(const Intrusive_set_column& c1, const Intrusive_set_column& c2){ - if (&c1 == &c2) return false; - - if constexpr (Master_matrix::Option_list::is_z2){ - return c1.column_ < c2.column_; - } else { - auto it1 = c1.column_.begin(); - auto it2 = c2.column_.begin(); - while (it1 != c1.column_.end() && it2 != c2.column_.end()) { - if (it1->get_row_index() != it2->get_row_index()) - return it1->get_row_index() < it2->get_row_index(); - if (it1->get_element() != it2->get_element()) - return it1->get_element() < it2->get_element(); - ++it1; ++it2; - } - return it2 != c2.column_.end(); - } - } - - // void set_operators(Field_operators* operators){ operators_ = operators; } - - //Disabled with row access. - Intrusive_set_column& operator=(const Intrusive_set_column& other); - - friend void swap(Intrusive_set_column& col1, Intrusive_set_column& col2){ - swap(static_cast(col1), - static_cast(col2)); - swap(static_cast(col1), - static_cast(col2)); - swap(static_cast(col1), - static_cast(col2)); - col1.column_.swap(col2.column_); - std::swap(col1.operators_, col2.operators_); - std::swap(col1.cellPool_, col2.cellPool_); - } - -private: - using ra_opt = typename Master_matrix::Row_access_option; - using dim_opt = typename Master_matrix::Column_dimension_option; - using chain_opt = typename Master_matrix::Chain_column_option; - - //Cloner object function for boost intrusive container - struct new_cloner - { - new_cloner(Cell_constructor* cellPool) : cellPool_(cellPool) - {}; - - Cell *operator()(const Cell& clone_this) - { - return cellPool_->construct(clone_this); - } - - Cell_constructor* cellPool_; - }; - - //The disposer object function for boost intrusive container - struct delete_disposer - { - delete_disposer(){}; - delete_disposer(Intrusive_set_column* col) : col_(col) - {}; - - void operator()(Cell *delete_this) - { - if constexpr (Master_matrix::Option_list::has_row_access) - col_->unlink(delete_this); - col_->cellPool_->destroy(delete_this); - } - - Intrusive_set_column* col_; - }; - - Column_type column_; - Field_operators* operators_; - Cell_constructor* cellPool_; - - void _delete_cell(iterator& it); - void _insert_cell(const Field_element_type& value, id_index rowIndex, const iterator& position); - void _insert_cell(id_index rowIndex, const iterator& position); - template - bool _add(const Cell_range& column); - template - bool _multiply_and_add(const Field_element_type& val, const Cell_range& column); - template - bool _multiply_and_add(const Cell_range& column, const Field_element_type& val); - - void _verifyCellConstructor(){ - if (cellPool_ == nullptr){ - if constexpr (std::is_same_v >){ - cellPool_ = &Master_matrix::defaultCellConstructor; - } else { - throw std::invalid_argument("Cell constructor pointer cannot be null."); - } - } - } +/** + * @brief Column class following the @ref PersistenceMatrixColumn concept. + * + * Column based on a intrusive set structure. The cells are ordered by row index and only non-zero values + * are stored uniquely in the underlying container. + * + * @tparam Master_matrix An instanciation of @ref Matrix from which all types and options are deduced. + * @tparam Cell_constructor Factory of @ref Cell classes. + */ +template > +class Intrusive_set_column : public Master_matrix::Row_access_option, + public Master_matrix::Column_dimension_option, + public Master_matrix::Chain_column_option +{ + public: + using Master = Master_matrix; + using Field_operators = typename Master_matrix::Field_operators; + using Field_element_type = typename Master_matrix::element_type; + using index = typename Master_matrix::index; + using id_index = typename Master_matrix::id_index; + using dimension_type = typename Master_matrix::dimension_type; + + using Cell = typename Master_matrix::Cell_type; + using Column_type = + boost::intrusive::set, + boost::intrusive::base_hook >; + using iterator = typename Column_type::iterator; + using const_iterator = typename Column_type::const_iterator; + using reverse_iterator = typename Column_type::reverse_iterator; + using const_reverse_iterator = typename Column_type::const_reverse_iterator; + + Intrusive_set_column(Field_operators* operators = nullptr, Cell_constructor* cellConstructor = nullptr); + template + Intrusive_set_column(const Container_type& nonZeroRowIndices, + Field_operators* operators, + Cell_constructor* cellConstructor); + template + Intrusive_set_column(index columnIndex, + const Container_type& nonZeroRowIndices, + Row_container_type* rowContainer, + Field_operators* operators, + Cell_constructor* cellConstructor); + template + Intrusive_set_column(const Container_type& nonZeroChainRowIndices, + dimension_type dimension, + Field_operators* operators, + Cell_constructor* cellConstructor); + template + Intrusive_set_column(index columnIndex, + const Container_type& nonZeroChainRowIndices, + dimension_type dimension, + Row_container_type* rowContainer, + Field_operators* operators, + Cell_constructor* cellConstructor); + Intrusive_set_column(const Intrusive_set_column& column, + Field_operators* operators = nullptr, + Cell_constructor* cellConstructor = nullptr); + template + Intrusive_set_column(const Intrusive_set_column& column, + index columnIndex, + Row_container_type* rowContainer, + Field_operators* operators = nullptr, + Cell_constructor* cellConstructor = nullptr); + Intrusive_set_column(Intrusive_set_column&& column) noexcept; + ~Intrusive_set_column(); + + std::vector get_content(int columnLength = -1) const; + bool is_non_zero(id_index rowIndex) const; + bool is_empty() const; + std::size_t size() const; + + template + void reorder(const Map_type& valueMap, [[maybe_unused]] index columnIndex = -1); + void clear(); + void clear(id_index rowIndex); + + id_index get_pivot() const; + Field_element_type get_pivot_value() const; + + iterator begin() noexcept; + const_iterator begin() const noexcept; + iterator end() noexcept; + const_iterator end() const noexcept; + reverse_iterator rbegin() noexcept; + const_reverse_iterator rbegin() const noexcept; + reverse_iterator rend() noexcept; + const_reverse_iterator rend() const noexcept; + + template + Intrusive_set_column& operator+=(const Cell_range& column); + Intrusive_set_column& operator+=(Intrusive_set_column& column); + + Intrusive_set_column& operator*=(unsigned int v); + + // this = v * this + column + template + Intrusive_set_column& multiply_and_add(const Field_element_type& val, const Cell_range& column); + Intrusive_set_column& multiply_and_add(const Field_element_type& val, Intrusive_set_column& column); + // this = this + column * v + template + Intrusive_set_column& multiply_and_add(const Cell_range& column, const Field_element_type& val); + Intrusive_set_column& multiply_and_add(Intrusive_set_column& column, const Field_element_type& val); + + friend bool operator==(const Intrusive_set_column& c1, const Intrusive_set_column& c2) { + if (&c1 == &c2) return true; + + if constexpr (Master_matrix::Option_list::is_z2) { + return c1.column_ == c2.column_; + } else { + auto it1 = c1.column_.begin(); + auto it2 = c2.column_.begin(); + if (c1.column_.size() != c2.column_.size()) return false; + while (it1 != c1.column_.end() && it2 != c2.column_.end()) { + if (it1->get_row_index() != it2->get_row_index() || it1->get_element() != it2->get_element()) return false; + ++it1; + ++it2; + } + return true; + } + } + friend bool operator<(const Intrusive_set_column& c1, const Intrusive_set_column& c2) { + if (&c1 == &c2) return false; + + if constexpr (Master_matrix::Option_list::is_z2) { + return c1.column_ < c2.column_; + } else { + auto it1 = c1.column_.begin(); + auto it2 = c2.column_.begin(); + while (it1 != c1.column_.end() && it2 != c2.column_.end()) { + if (it1->get_row_index() != it2->get_row_index()) return it1->get_row_index() < it2->get_row_index(); + if (it1->get_element() != it2->get_element()) return it1->get_element() < it2->get_element(); + ++it1; + ++it2; + } + return it2 != c2.column_.end(); + } + } + + // void set_operators(Field_operators* operators){ operators_ = operators; } + + // Disabled with row access. + Intrusive_set_column& operator=(const Intrusive_set_column& other); + + friend void swap(Intrusive_set_column& col1, Intrusive_set_column& col2) { + swap(static_cast(col1), + static_cast(col2)); + swap(static_cast(col1), + static_cast(col2)); + swap(static_cast(col1), + static_cast(col2)); + col1.column_.swap(col2.column_); + std::swap(col1.operators_, col2.operators_); + std::swap(col1.cellPool_, col2.cellPool_); + } + + private: + using ra_opt = typename Master_matrix::Row_access_option; + using dim_opt = typename Master_matrix::Column_dimension_option; + using chain_opt = typename Master_matrix::Chain_column_option; + + // Cloner object function for boost intrusive container + struct new_cloner { + new_cloner(Cell_constructor* cellPool) : cellPool_(cellPool){}; + + Cell* operator()(const Cell& clone_this) { return cellPool_->construct(clone_this); } + + Cell_constructor* cellPool_; + }; + + // The disposer object function for boost intrusive container + struct delete_disposer { + delete_disposer(){}; + delete_disposer(Intrusive_set_column* col) : col_(col){}; + + void operator()(Cell* delete_this) { + if constexpr (Master_matrix::Option_list::has_row_access) col_->unlink(delete_this); + col_->cellPool_->destroy(delete_this); + } + + Intrusive_set_column* col_; + }; + + Column_type column_; + Field_operators* operators_; + Cell_constructor* cellPool_; + + void _delete_cell(iterator& it); + void _insert_cell(const Field_element_type& value, id_index rowIndex, const iterator& position); + void _insert_cell(id_index rowIndex, const iterator& position); + template + bool _add(const Cell_range& column); + template + bool _multiply_and_add(const Field_element_type& val, const Cell_range& column); + template + bool _multiply_and_add(const Cell_range& column, const Field_element_type& val); + + void _verifyCellConstructor() { + if (cellPool_ == nullptr) { + if constexpr (std::is_same_v >) { + cellPool_ = &Master_matrix::defaultCellConstructor; + } else { + throw std::invalid_argument("Cell constructor pointer cannot be null."); + } + } + } }; -template -inline Intrusive_set_column::Intrusive_set_column(Field_operators* operators, Cell_constructor* cellConstructor) : ra_opt(), dim_opt(), chain_opt(), operators_(operators), - cellPool_(cellConstructor) -{ - if (operators_ == nullptr && cellPool_ == nullptr) return; //to allow default constructor which gives a dummy column - _verifyCellConstructor(); -} - -template -template -inline Intrusive_set_column::Intrusive_set_column(const Container_type &nonZeroRowIndices, Field_operators* operators, Cell_constructor* cellConstructor) - : ra_opt(), - dim_opt(nonZeroRowIndices.size() == 0 ? 0 : nonZeroRowIndices.size() - 1), - chain_opt(), - operators_(operators), - cellPool_(cellConstructor) -{ - static_assert(!Master_matrix::isNonBasic || Master_matrix::Option_list::is_of_boundary_type, - "Constructor not available for chain columns, please specify the dimension of the chain."); - - _verifyCellConstructor(); - - if constexpr (Master_matrix::Option_list::is_z2){ - for (id_index id : nonZeroRowIndices){ - _insert_cell(id, column_.end()); - } - } else { - for (const auto& p : nonZeroRowIndices){ - _insert_cell(operators_->get_value(p.second), p.first, column_.end()); - } - } -} - -template -template -inline Intrusive_set_column::Intrusive_set_column( - index columnIndex, const Container_type &nonZeroRowIndices, Row_container_type *rowContainer, Field_operators* operators, Cell_constructor* cellConstructor) - : ra_opt(columnIndex, rowContainer), - dim_opt(nonZeroRowIndices.size() == 0 ? 0 : nonZeroRowIndices.size() - 1), - chain_opt([&]{ - if constexpr (Master_matrix::Option_list::is_z2){ - return nonZeroRowIndices.begin() == nonZeroRowIndices.end() ? -1 : *std::prev(nonZeroRowIndices.end()); - } else { - return nonZeroRowIndices.begin() == nonZeroRowIndices.end() ? -1 : std::prev(nonZeroRowIndices.end())->first; - } - }()), - operators_(operators), - cellPool_(cellConstructor) -{ - static_assert(!Master_matrix::isNonBasic || Master_matrix::Option_list::is_of_boundary_type, - "Constructor not available for chain columns, please specify the dimension of the chain."); - - _verifyCellConstructor(); - - if constexpr (Master_matrix::Option_list::is_z2){ - for (id_index id : nonZeroRowIndices){ - _insert_cell(id, column_.end()); - } - } else { - for (const auto& p : nonZeroRowIndices){ - _insert_cell(operators_->get_value(p.second), p.first, column_.end()); - } - } -} - -template -template -inline Intrusive_set_column::Intrusive_set_column( - const Container_type &nonZeroRowIndices, dimension_type dimension, Field_operators* operators, Cell_constructor* cellConstructor) - : ra_opt(), - dim_opt(dimension), - chain_opt([&]{ - if constexpr (Master_matrix::Option_list::is_z2){ - return nonZeroRowIndices.begin() == nonZeroRowIndices.end() ? -1 : *std::prev(nonZeroRowIndices.end()); - } else { - return nonZeroRowIndices.begin() == nonZeroRowIndices.end() ? -1 : std::prev(nonZeroRowIndices.end())->first; - } - }()), - operators_(operators), - cellPool_(cellConstructor) -{ - _verifyCellConstructor(); - - if constexpr (Master_matrix::Option_list::is_z2){ - for (id_index id : nonZeroRowIndices){ - _insert_cell(id, column_.end()); - } - } else { - for (const auto& p : nonZeroRowIndices){ - _insert_cell(operators_->get_value(p.second), p.first, column_.end()); - } - } -} - -template -template -inline Intrusive_set_column::Intrusive_set_column( - index columnIndex, const Container_type &nonZeroRowIndices, dimension_type dimension, Row_container_type *rowContainer, Field_operators* operators, Cell_constructor* cellConstructor) - : ra_opt(columnIndex, rowContainer), - dim_opt(dimension), - chain_opt([&]{ - if constexpr (Master_matrix::Option_list::is_z2){ - return nonZeroRowIndices.begin() == nonZeroRowIndices.end() ? -1 : *std::prev(nonZeroRowIndices.end()); - } else { - return nonZeroRowIndices.begin() == nonZeroRowIndices.end() ? -1 : std::prev(nonZeroRowIndices.end())->first; - } - }()), - operators_(operators), - cellPool_(cellConstructor) -{ - _verifyCellConstructor(); - - if constexpr (Master_matrix::Option_list::is_z2){ - for (id_index id : nonZeroRowIndices){ - _insert_cell(id, column_.end()); - } - } else { - for (const auto& p : nonZeroRowIndices){ - _insert_cell(operators_->get_value(p.second), p.first, column_.end()); - } - } -} - -template -inline Intrusive_set_column::Intrusive_set_column(const Intrusive_set_column &column, Field_operators* operators, Cell_constructor* cellConstructor) - : ra_opt(), - dim_opt(static_cast(column)), - chain_opt(static_cast(column)), - operators_(operators == nullptr ? column.operators_ : operators), - cellPool_(cellConstructor == nullptr ? column.cellPool_ : cellConstructor) -{ - static_assert(!Master_matrix::Option_list::has_row_access, - "Simple copy constructor not available when row access option enabled. Please specify the new column index and the row container."); - - column_.clone_from(column.column_, new_cloner(cellPool_), delete_disposer(this)); -} - -template -template -inline Intrusive_set_column::Intrusive_set_column( - const Intrusive_set_column &column, index columnIndex, Row_container_type *rowContainer, Field_operators* operators, Cell_constructor* cellConstructor) - : ra_opt(columnIndex, rowContainer), - dim_opt(static_cast(column)), - chain_opt(static_cast(column)), - operators_(operators == nullptr ? column.operators_ : operators), - cellPool_(cellConstructor == nullptr ? column.cellPool_ : cellConstructor) -{ - for (const Cell& cell : column.column_){ - if constexpr (Master_matrix::Option_list::is_z2){ - _insert_cell(cell.get_row_index(), column_.end()); - } else { - _insert_cell(cell.get_element(), cell.get_row_index(), column_.end()); - } - } -} - -template -inline Intrusive_set_column::Intrusive_set_column(Intrusive_set_column &&column) noexcept - : ra_opt(std::move(static_cast(column))), - dim_opt(std::move(static_cast(column))), - chain_opt(std::move(static_cast(column))), - column_(std::move(column.column_)), - operators_(std::exchange(column.operators_, nullptr)), - cellPool_(std::exchange(column.cellPool_, nullptr)) +template +inline Intrusive_set_column::Intrusive_set_column(Field_operators* operators, + Cell_constructor* cellConstructor) + : ra_opt(), dim_opt(), chain_opt(), operators_(operators), cellPool_(cellConstructor) +{ + if (operators_ == nullptr && cellPool_ == nullptr) return; //to allow default constructor which gives a dummy column + _verifyCellConstructor(); +} + +template +template +inline Intrusive_set_column::Intrusive_set_column( + const Container_type& nonZeroRowIndices, Field_operators* operators, Cell_constructor* cellConstructor) + : ra_opt(), + dim_opt(nonZeroRowIndices.size() == 0 ? 0 : nonZeroRowIndices.size() - 1), + chain_opt(), + operators_(operators), + cellPool_(cellConstructor) +{ + static_assert(!Master_matrix::isNonBasic || Master_matrix::Option_list::is_of_boundary_type, + "Constructor not available for chain columns, please specify the dimension of the chain."); + + _verifyCellConstructor(); + + if constexpr (Master_matrix::Option_list::is_z2) { + for (id_index id : nonZeroRowIndices) { + _insert_cell(id, column_.end()); + } + } else { + for (const auto& p : nonZeroRowIndices) { + _insert_cell(operators_->get_value(p.second), p.first, column_.end()); + } + } +} + +template +template +inline Intrusive_set_column::Intrusive_set_column( + index columnIndex, + const Container_type& nonZeroRowIndices, + Row_container_type* rowContainer, + Field_operators* operators, + Cell_constructor* cellConstructor) + : ra_opt(columnIndex, rowContainer), + dim_opt(nonZeroRowIndices.size() == 0 ? 0 : nonZeroRowIndices.size() - 1), + chain_opt([&] { + if constexpr (Master_matrix::Option_list::is_z2) { + return nonZeroRowIndices.begin() == nonZeroRowIndices.end() ? -1 : *std::prev(nonZeroRowIndices.end()); + } else { + return nonZeroRowIndices.begin() == nonZeroRowIndices.end() ? -1 : std::prev(nonZeroRowIndices.end())->first; + } + }()), + operators_(operators), + cellPool_(cellConstructor) +{ + static_assert(!Master_matrix::isNonBasic || Master_matrix::Option_list::is_of_boundary_type, + "Constructor not available for chain columns, please specify the dimension of the chain."); + + _verifyCellConstructor(); + + if constexpr (Master_matrix::Option_list::is_z2) { + for (id_index id : nonZeroRowIndices) { + _insert_cell(id, column_.end()); + } + } else { + for (const auto& p : nonZeroRowIndices) { + _insert_cell(operators_->get_value(p.second), p.first, column_.end()); + } + } +} + +template +template +inline Intrusive_set_column::Intrusive_set_column( + const Container_type& nonZeroRowIndices, + dimension_type dimension, + Field_operators* operators, + Cell_constructor* cellConstructor) + : ra_opt(), + dim_opt(dimension), + chain_opt([&] { + if constexpr (Master_matrix::Option_list::is_z2) { + return nonZeroRowIndices.begin() == nonZeroRowIndices.end() ? -1 : *std::prev(nonZeroRowIndices.end()); + } else { + return nonZeroRowIndices.begin() == nonZeroRowIndices.end() ? -1 : std::prev(nonZeroRowIndices.end())->first; + } + }()), + operators_(operators), + cellPool_(cellConstructor) +{ + _verifyCellConstructor(); + + if constexpr (Master_matrix::Option_list::is_z2) { + for (id_index id : nonZeroRowIndices) { + _insert_cell(id, column_.end()); + } + } else { + for (const auto& p : nonZeroRowIndices) { + _insert_cell(operators_->get_value(p.second), p.first, column_.end()); + } + } +} + +template +template +inline Intrusive_set_column::Intrusive_set_column( + index columnIndex, + const Container_type& nonZeroRowIndices, + dimension_type dimension, + Row_container_type* rowContainer, + Field_operators* operators, + Cell_constructor* cellConstructor) + : ra_opt(columnIndex, rowContainer), + dim_opt(dimension), + chain_opt([&] { + if constexpr (Master_matrix::Option_list::is_z2) { + return nonZeroRowIndices.begin() == nonZeroRowIndices.end() ? -1 : *std::prev(nonZeroRowIndices.end()); + } else { + return nonZeroRowIndices.begin() == nonZeroRowIndices.end() ? -1 : std::prev(nonZeroRowIndices.end())->first; + } + }()), + operators_(operators), + cellPool_(cellConstructor) { + _verifyCellConstructor(); + + if constexpr (Master_matrix::Option_list::is_z2) { + for (id_index id : nonZeroRowIndices) { + _insert_cell(id, column_.end()); + } + } else { + for (const auto& p : nonZeroRowIndices) { + _insert_cell(operators_->get_value(p.second), p.first, column_.end()); + } + } +} + +template +inline Intrusive_set_column::Intrusive_set_column(const Intrusive_set_column& column, + Field_operators* operators, + Cell_constructor* cellConstructor) + : ra_opt(), + dim_opt(static_cast(column)), + chain_opt(static_cast(column)), + operators_(operators == nullptr ? column.operators_ : operators), + cellPool_(cellConstructor == nullptr ? column.cellPool_ : cellConstructor) +{ + static_assert(!Master_matrix::Option_list::has_row_access, + "Simple copy constructor not available when row access option enabled. Please specify the new column " + "index and the row container."); + + column_.clone_from(column.column_, new_cloner(cellPool_), delete_disposer(this)); +} + +template +template +inline Intrusive_set_column::Intrusive_set_column(const Intrusive_set_column& column, + index columnIndex, + Row_container_type* rowContainer, + Field_operators* operators, + Cell_constructor* cellConstructor) + : ra_opt(columnIndex, rowContainer), + dim_opt(static_cast(column)), + chain_opt(static_cast(column)), + operators_(operators == nullptr ? column.operators_ : operators), + cellPool_(cellConstructor == nullptr ? column.cellPool_ : cellConstructor) +{ + for (const Cell& cell : column.column_) { + if constexpr (Master_matrix::Option_list::is_z2) { + _insert_cell(cell.get_row_index(), column_.end()); + } else { + _insert_cell(cell.get_element(), cell.get_row_index(), column_.end()); + } + } +} + +template +inline Intrusive_set_column::Intrusive_set_column( + Intrusive_set_column&& column) noexcept + : ra_opt(std::move(static_cast(column))), + dim_opt(std::move(static_cast(column))), + chain_opt(std::move(static_cast(column))), + column_(std::move(column.column_)), + operators_(std::exchange(column.operators_, nullptr)), + cellPool_(std::exchange(column.cellPool_, nullptr)) {} -template -inline Intrusive_set_column::~Intrusive_set_column() +template +inline Intrusive_set_column::~Intrusive_set_column() { - column_.clear_and_dispose(delete_disposer(this)); + column_.clear_and_dispose(delete_disposer(this)); } -template -inline std::vector::Field_element_type> -Intrusive_set_column::get_content(int columnLength) const +template +inline std::vector::Field_element_type> +Intrusive_set_column::get_content(int columnLength) const { - if (columnLength < 0 && column_.size() > 0) columnLength = column_.rbegin()->get_row_index() + 1; - else if (columnLength < 0) return std::vector(); + if (columnLength < 0 && column_.size() > 0) + columnLength = column_.rbegin()->get_row_index() + 1; + else if (columnLength < 0) + return std::vector(); - std::vector container(columnLength); - for (auto it = column_.begin(); it != column_.end() && it->get_row_index() < static_cast(columnLength); ++it){ - if constexpr (Master_matrix::Option_list::is_z2){ - container[it->get_row_index()] = 1; - } else { - container[it->get_row_index()] = it->get_element(); - } - } - return container; + std::vector container(columnLength); + for (auto it = column_.begin(); it != column_.end() && it->get_row_index() < static_cast(columnLength); + ++it) { + if constexpr (Master_matrix::Option_list::is_z2) { + container[it->get_row_index()] = 1; + } else { + container[it->get_row_index()] = it->get_element(); + } + } + return container; } -template -inline bool Intrusive_set_column::is_non_zero(id_index rowIndex) const +template +inline bool Intrusive_set_column::is_non_zero(id_index rowIndex) const { - return column_.find(Cell(rowIndex)) != column_.end(); + return column_.find(Cell(rowIndex)) != column_.end(); } -template -inline bool Intrusive_set_column::is_empty() const +template +inline bool Intrusive_set_column::is_empty() const { - return column_.empty(); + return column_.empty(); } -template -inline std::size_t Intrusive_set_column::size() const{ - return column_.size(); +template +inline std::size_t Intrusive_set_column::size() const +{ + return column_.size(); } -template -template -inline void Intrusive_set_column::reorder(const Map_type &valueMap, [[maybe_unused]] index columnIndex) +template +template +inline void Intrusive_set_column::reorder(const Map_type& valueMap, + [[maybe_unused]] index columnIndex) { - static_assert(!Master_matrix::isNonBasic || Master_matrix::Option_list::is_of_boundary_type, - "Method not available for chain columns."); + static_assert(!Master_matrix::isNonBasic || Master_matrix::Option_list::is_of_boundary_type, + "Method not available for chain columns."); - Column_type newSet; + Column_type newSet; - if constexpr (Master_matrix::Option_list::has_row_access) { - for (auto it = column_.begin(); it != column_.end(); ) { - Cell *new_cell = cellPool_->construct(columnIndex == static_cast(-1) ? ra_opt::columnIndex_ : columnIndex, valueMap.at(it->get_row_index())); - if constexpr (!Master_matrix::Option_list::is_z2){ - new_cell->set_element(it->get_element()); - } - newSet.insert(newSet.end(), *new_cell); - _delete_cell(it); //increases it - if constexpr (Master_matrix::Option_list::has_intrusive_rows) //intrusive list - ra_opt::insert_cell(new_cell->get_row_index(), new_cell); - } + if constexpr (Master_matrix::Option_list::has_row_access) { + for (auto it = column_.begin(); it != column_.end();) { + Cell* new_cell = cellPool_->construct(columnIndex == static_cast(-1) ? ra_opt::columnIndex_ : columnIndex, + valueMap.at(it->get_row_index())); + if constexpr (!Master_matrix::Option_list::is_z2) { + new_cell->set_element(it->get_element()); + } + newSet.insert(newSet.end(), *new_cell); + _delete_cell(it); // increases it + if constexpr (Master_matrix::Option_list::has_intrusive_rows) // intrusive list + ra_opt::insert_cell(new_cell->get_row_index(), new_cell); + } - //when row is a set, all cells have to be deleted first, to avoid colliding when inserting - if constexpr (!Master_matrix::Option_list::has_intrusive_rows){ //set - for (Cell& cell : newSet) { - ra_opt::insert_cell(cell.get_row_index(), &cell); - } - } - } else { - for (auto it = column_.begin(); it != column_.end(); ) { - Cell *new_cell = cellPool_->construct(valueMap.at(it->get_row_index())); - if constexpr (!Master_matrix::Option_list::is_z2){ - new_cell->set_element(it->get_element()); - } - newSet.insert(newSet.end(), *new_cell); - _delete_cell(it); //increases it - } - } + // when row is a set, all cells have to be deleted first, to avoid colliding when inserting + if constexpr (!Master_matrix::Option_list::has_intrusive_rows) { // set + for (Cell& cell : newSet) { + ra_opt::insert_cell(cell.get_row_index(), &cell); + } + } + } else { + for (auto it = column_.begin(); it != column_.end();) { + Cell* new_cell = cellPool_->construct(valueMap.at(it->get_row_index())); + if constexpr (!Master_matrix::Option_list::is_z2) { + new_cell->set_element(it->get_element()); + } + newSet.insert(newSet.end(), *new_cell); + _delete_cell(it); // increases it + } + } - column_.swap(newSet); + column_.swap(newSet); } -template -inline void Intrusive_set_column::clear() +template +inline void Intrusive_set_column::clear() { - static_assert(!Master_matrix::isNonBasic || Master_matrix::Option_list::is_of_boundary_type, - "Method not available for chain columns as a base element should not be empty."); + static_assert(!Master_matrix::isNonBasic || Master_matrix::Option_list::is_of_boundary_type, + "Method not available for chain columns as a base element should not be empty."); - column_.clear_and_dispose(delete_disposer(this)); + column_.clear_and_dispose(delete_disposer(this)); } -template -inline void Intrusive_set_column::clear(id_index rowIndex) +template +inline void Intrusive_set_column::clear(id_index rowIndex) { - static_assert(!Master_matrix::isNonBasic || Master_matrix::Option_list::is_of_boundary_type, - "Method not available for chain columns."); + static_assert(!Master_matrix::isNonBasic || Master_matrix::Option_list::is_of_boundary_type, + "Method not available for chain columns."); - auto it = column_.find(Cell(rowIndex)); - if (it != column_.end()){ - _delete_cell(it); - } + auto it = column_.find(Cell(rowIndex)); + if (it != column_.end()) { + _delete_cell(it); + } } -template -inline typename Intrusive_set_column::id_index Intrusive_set_column::get_pivot() const +template +inline typename Intrusive_set_column::id_index +Intrusive_set_column::get_pivot() const { - static_assert(Master_matrix::isNonBasic, "Method not available for base columns."); //could technically be, but is the notion usefull then? + static_assert(Master_matrix::isNonBasic, + "Method not available for base columns."); // could technically be, but is the notion usefull then? - if constexpr (Master_matrix::Option_list::is_of_boundary_type){ - if (column_.empty()) return -1; - return column_.rbegin()->get_row_index(); - } else { - return chain_opt::get_pivot(); - } + if constexpr (Master_matrix::Option_list::is_of_boundary_type) { + if (column_.empty()) return -1; + return column_.rbegin()->get_row_index(); + } else { + return chain_opt::get_pivot(); + } } -template -inline typename Intrusive_set_column::Field_element_type Intrusive_set_column::get_pivot_value() const +template +inline typename Intrusive_set_column::Field_element_type +Intrusive_set_column::get_pivot_value() const { - static_assert(Master_matrix::isNonBasic, "Method not available for base columns."); //could technically be, but is the notion usefull then? + static_assert(Master_matrix::isNonBasic, + "Method not available for base columns."); // could technically be, but is the notion usefull then? - if constexpr (Master_matrix::Option_list::is_z2){ - return 1; - } else { - if constexpr (Master_matrix::Option_list::is_of_boundary_type){ - if (column_.empty()) return 0; - return column_.rbegin()->get_element(); - } else { - if (chain_opt::get_pivot() == -1) return 0; - auto it = column_.find(Cell(chain_opt::get_pivot())); - GUDHI_CHECK(it != column_.end(), "Pivot not found only if the column was misused."); - return it->get_element(); - } - } + if constexpr (Master_matrix::Option_list::is_z2) { + return 1; + } else { + if constexpr (Master_matrix::Option_list::is_of_boundary_type) { + if (column_.empty()) return 0; + return column_.rbegin()->get_element(); + } else { + if (chain_opt::get_pivot() == -1) return 0; + auto it = column_.find(Cell(chain_opt::get_pivot())); + GUDHI_CHECK(it != column_.end(), "Pivot not found only if the column was misused."); + return it->get_element(); + } + } } -template -inline typename Intrusive_set_column::iterator -Intrusive_set_column::begin() noexcept +template +inline typename Intrusive_set_column::iterator +Intrusive_set_column::begin() noexcept { - return column_.begin(); + return column_.begin(); } -template -inline typename Intrusive_set_column::const_iterator -Intrusive_set_column::begin() const noexcept +template +inline typename Intrusive_set_column::const_iterator +Intrusive_set_column::begin() const noexcept { - return column_.begin(); + return column_.begin(); } -template -inline typename Intrusive_set_column::iterator -Intrusive_set_column::end() noexcept +template +inline typename Intrusive_set_column::iterator +Intrusive_set_column::end() noexcept { - return column_.end(); + return column_.end(); } -template -inline typename Intrusive_set_column::const_iterator -Intrusive_set_column::end() const noexcept +template +inline typename Intrusive_set_column::const_iterator +Intrusive_set_column::end() const noexcept { - return column_.end(); + return column_.end(); } -template -inline typename Intrusive_set_column::reverse_iterator -Intrusive_set_column::rbegin() noexcept +template +inline typename Intrusive_set_column::reverse_iterator +Intrusive_set_column::rbegin() noexcept { - return column_.rbegin(); + return column_.rbegin(); } -template -inline typename Intrusive_set_column::const_reverse_iterator -Intrusive_set_column::rbegin() const noexcept +template +inline typename Intrusive_set_column::const_reverse_iterator +Intrusive_set_column::rbegin() const noexcept { - return column_.rbegin(); + return column_.rbegin(); } -template -inline typename Intrusive_set_column::reverse_iterator -Intrusive_set_column::rend() noexcept +template +inline typename Intrusive_set_column::reverse_iterator +Intrusive_set_column::rend() noexcept { - return column_.rend(); + return column_.rend(); } -template -inline typename Intrusive_set_column::const_reverse_iterator -Intrusive_set_column::rend() const noexcept +template +inline typename Intrusive_set_column::const_reverse_iterator +Intrusive_set_column::rend() const noexcept { - return column_.rend(); + return column_.rend(); } -template -template -inline Intrusive_set_column & -Intrusive_set_column::operator+=(const Cell_range &column) +template +template +inline Intrusive_set_column& +Intrusive_set_column::operator+=(const Cell_range& column) { - static_assert((!Master_matrix::isNonBasic || std::is_same_v), - "For boundary columns, the range has to be a column of same type to help ensure the validity of the base element."); //could be removed, if we give the responsability to the user. - static_assert((!Master_matrix::isNonBasic || Master_matrix::Option_list::is_of_boundary_type), - "For chain columns, the given column cannot be constant."); + static_assert((!Master_matrix::isNonBasic || std::is_same_v), + "For boundary columns, the range has to be a column of same type to help ensure the validity of the " + "base element."); // could be removed, if we give the responsability to the user. + static_assert((!Master_matrix::isNonBasic || Master_matrix::Option_list::is_of_boundary_type), + "For chain columns, the given column cannot be constant."); - _add(column); + _add(column); - return *this; + return *this; } -template -inline Intrusive_set_column & -Intrusive_set_column::operator+=(Intrusive_set_column &column) +template +inline Intrusive_set_column& +Intrusive_set_column::operator+=(Intrusive_set_column& column) { - if constexpr (Master_matrix::isNonBasic && !Master_matrix::Option_list::is_of_boundary_type){ - //assumes that the addition never zeros out this column. - if (_add(column)) { - chain_opt::swap_pivots(column); - dim_opt::swap_dimension(column); - } - } else { - _add(column); - } - - return *this; + if constexpr (Master_matrix::isNonBasic && !Master_matrix::Option_list::is_of_boundary_type) { + // assumes that the addition never zeros out this column. + if (_add(column)) { + chain_opt::swap_pivots(column); + dim_opt::swap_dimension(column); + } + } else { + _add(column); + } + + return *this; } -template -inline Intrusive_set_column & -Intrusive_set_column::operator*=(unsigned int v) +template +inline Intrusive_set_column& +Intrusive_set_column::operator*=(unsigned int v) { - if constexpr (Master_matrix::Option_list::is_z2){ - if (v % 2 == 0){ - if constexpr (Master_matrix::isNonBasic && !Master_matrix::Option_list::is_of_boundary_type){ - throw std::invalid_argument("A chain column should not be multiplied by 0."); - } else { - clear(); - } - } - } else { - Field_element_type val = operators_->get_value(v); + if constexpr (Master_matrix::Option_list::is_z2) { + if (v % 2 == 0) { + if constexpr (Master_matrix::isNonBasic && !Master_matrix::Option_list::is_of_boundary_type) { + throw std::invalid_argument("A chain column should not be multiplied by 0."); + } else { + clear(); + } + } + } else { + Field_element_type val = operators_->get_value(v); - if (val == Field_operators::get_additive_identity()) { - if constexpr (Master_matrix::isNonBasic && !Master_matrix::Option_list::is_of_boundary_type){ - throw std::invalid_argument("A chain column should not be multiplied by 0."); - } else { - clear(); - } - return *this; - } + if (val == Field_operators::get_additive_identity()) { + if constexpr (Master_matrix::isNonBasic && !Master_matrix::Option_list::is_of_boundary_type) { + throw std::invalid_argument("A chain column should not be multiplied by 0."); + } else { + clear(); + } + return *this; + } - if (val == Field_operators::get_multiplicative_identity()) return *this; - - for (Cell& cell : column_){ - cell.get_element() = operators_->multiply(cell.get_element(), val); - if constexpr (Master_matrix::Option_list::has_row_access) - ra_opt::update_cell(cell); - } - } - - return *this; -} - -template -template -inline Intrusive_set_column & -Intrusive_set_column::multiply_and_add(const Field_element_type& val, const Cell_range& column) -{ - static_assert((!Master_matrix::isNonBasic || std::is_same_v), - "For boundary columns, the range has to be a column of same type to help ensure the validity of the base element."); //could be removed, if we give the responsability to the user. - static_assert((!Master_matrix::isNonBasic || Master_matrix::Option_list::is_of_boundary_type), - "For chain columns, the given column cannot be constant."); - - if constexpr (Master_matrix::Option_list::is_z2){ - if (val){ - _add(column); - } else { - clear(); - _add(column); - } - } else { - _multiply_and_add(val, column); - } - - return *this; -} - -template -inline Intrusive_set_column & -Intrusive_set_column::multiply_and_add(const Field_element_type& val, Intrusive_set_column& column) -{ - if constexpr (Master_matrix::isNonBasic && !Master_matrix::Option_list::is_of_boundary_type){ - //assumes that the addition never zeros out this column. - if constexpr (Master_matrix::Option_list::is_z2){ - if (val){ - if (_add(column)){ - chain_opt::swap_pivots(column); - dim_opt::swap_dimension(column); - } - } else { - throw std::invalid_argument("A chain column should not be multiplied by 0."); - } - } else { - if (_multiply_and_add(val, column)){ - chain_opt::swap_pivots(column); - dim_opt::swap_dimension(column); - } - } - } else { - if constexpr (Master_matrix::Option_list::is_z2){ - if (val){ - _add(column); - } else { - clear(); - _add(column); - } - } else { - _multiply_and_add(val, column); - } - } - - return *this; -} - -template -template -inline Intrusive_set_column & -Intrusive_set_column::multiply_and_add(const Cell_range& column, const Field_element_type& val) -{ - static_assert((!Master_matrix::isNonBasic || std::is_same_v), - "For boundary columns, the range has to be a column of same type to help ensure the validity of the base element."); //could be removed, if we give the responsability to the user. - static_assert((!Master_matrix::isNonBasic || Master_matrix::Option_list::is_of_boundary_type), - "For chain columns, the given column cannot be constant."); - - if constexpr (Master_matrix::Option_list::is_z2){ - if (val){ - _add(column); - } - } else { - _multiply_and_add(column, val); - } - - return *this; -} - -template -inline Intrusive_set_column & -Intrusive_set_column::multiply_and_add(Intrusive_set_column& column, const Field_element_type& val) -{ - if constexpr (Master_matrix::isNonBasic && !Master_matrix::Option_list::is_of_boundary_type){ - //assumes that the addition never zeros out this column. - if constexpr (Master_matrix::Option_list::is_z2){ - if (val){ - if (_add(column)){ - chain_opt::swap_pivots(column); - dim_opt::swap_dimension(column); - } - } - } else { - if (_multiply_and_add(column, val)){ - chain_opt::swap_pivots(column); - dim_opt::swap_dimension(column); - } - } - } else { - if constexpr (Master_matrix::Option_list::is_z2){ - if (val){ - _add(column); - } - } else { - _multiply_and_add(column, val); - } - } - - return *this; -} - -template -inline Intrusive_set_column & -Intrusive_set_column::operator=(const Intrusive_set_column& other) -{ - static_assert(!Master_matrix::Option_list::has_row_access, "= assignement not enabled with row access option."); - - dim_opt::operator=(other); - chain_opt::operator=(other); - - //order is important - column_.clear_and_dispose(delete_disposer(this)); - operators_ = other.operators_; - cellPool_ = other.cellPool_; - column_.clone_from(other.column_, new_cloner(cellPool_), delete_disposer(this)); - - return *this; -} - -template -inline void Intrusive_set_column::_delete_cell(iterator &it) -{ - it = column_.erase_and_dispose(it, delete_disposer(this)); -} - -template -inline void Intrusive_set_column::_insert_cell( - const Field_element_type &value, id_index rowIndex, const iterator &position) -{ - if constexpr (Master_matrix::Option_list::has_row_access){ - Cell *new_cell = cellPool_->construct(ra_opt::columnIndex_, rowIndex); - new_cell->set_element(value); - column_.insert(position, *new_cell); - ra_opt::insert_cell(rowIndex, new_cell); - } else { - Cell *new_cell = cellPool_->construct(rowIndex); - new_cell->set_element(value); - column_.insert(position, *new_cell); - } -} - -template -inline void Intrusive_set_column::_insert_cell( - id_index rowIndex, const iterator &position) -{ - if constexpr (Master_matrix::Option_list::has_row_access){ - Cell *new_cell = cellPool_->construct(ra_opt::columnIndex_, rowIndex); - column_.insert(position, *new_cell); - ra_opt::insert_cell(rowIndex, new_cell); - } else { - Cell *new_cell = cellPool_->construct(rowIndex); - column_.insert(position, *new_cell); - } -} - -template -template -inline bool Intrusive_set_column::_add(const Cell_range &column) -{ - bool pivotIsZeroed = false; - - for (const Cell &cell : column) { - auto it1 = column_.find(cell); - if (it1 != column_.end()) { - if constexpr (Master_matrix::Option_list::is_z2){ - if constexpr (Master_matrix::isNonBasic && !Master_matrix::Option_list::is_of_boundary_type){ - if (it1->get_row_index() == chain_opt::get_pivot()) pivotIsZeroed = true; - } - _delete_cell(it1); - } else { - it1->get_element() = operators_->add(it1->get_element(), cell.get_element()); - if (it1->get_element() == Field_operators::get_additive_identity()){ - if constexpr (Master_matrix::isNonBasic && !Master_matrix::Option_list::is_of_boundary_type){ - if (it1->get_row_index() == chain_opt::get_pivot()) pivotIsZeroed = true; - } - _delete_cell(it1); - } else { - if constexpr (Master_matrix::Option_list::has_row_access) - ra_opt::update_cell(*it1); - } - } - } else { - if constexpr (Master_matrix::Option_list::is_z2){ - _insert_cell(cell.get_row_index(), column_.end()); - } else { - _insert_cell(cell.get_element(), cell.get_row_index(), column_.end()); - } - } - } - - return pivotIsZeroed; -} - -template -template -inline bool Intrusive_set_column::_multiply_and_add(const Field_element_type& val, const Cell_range& column) -{ - bool pivotIsZeroed = false; - - if (val == 0u) { - if constexpr (Master_matrix::isNonBasic && !Master_matrix::Option_list::is_of_boundary_type){ - throw std::invalid_argument("A chain column should not be multiplied by 0."); - //this would not only mess up the base, but also the pivots stored. - } else { - clear(); - } - } - - //cannot use usual addition method, because all values of column_ have to be multiplied - - auto itTarget = column_.begin(); - auto itSource = column.begin(); - while (itTarget != column_.end() && itSource != column.end()) - { - if (itTarget->get_row_index() < itSource->get_row_index()) { - itTarget->get_element() = operators_->multiply(itTarget->get_element(), val); - if constexpr (Master_matrix::Option_list::has_row_access) - ra_opt::update_cell(*itTarget); - ++itTarget; - } else if (itTarget->get_row_index() > itSource->get_row_index()) { - _insert_cell(itSource->get_element(), itSource->get_row_index(), itTarget); - ++itSource; - } else { - itTarget->get_element() = operators_->multiply_and_add(itTarget->get_element(), val, itSource->get_element()); - if (itTarget->get_element() == Field_operators::get_additive_identity()){ - if constexpr (Master_matrix::isNonBasic && !Master_matrix::Option_list::is_of_boundary_type){ - if (itTarget->get_row_index() == chain_opt::get_pivot()) pivotIsZeroed = true; - } - _delete_cell(itTarget); - } else { - if constexpr (Master_matrix::Option_list::has_row_access) - ra_opt::update_cell(*itTarget); - ++itTarget; - } - ++itSource; - } - } - - while (itTarget != column_.end()){ - itTarget->get_element() = operators_->multiply(itTarget->get_element(), val); - if constexpr (Master_matrix::Option_list::has_row_access) - ra_opt::update_cell(*itTarget); - itTarget++; - } - - while (itSource != column.end()) { - _insert_cell(itSource->get_element(), itSource->get_row_index(), column_.end()); - ++itSource; - } - - return pivotIsZeroed; -} - -template -template -inline bool Intrusive_set_column::_multiply_and_add(const Cell_range& column, const Field_element_type& val) -{ - if (val == 0u) { - return false; - } - - bool pivotIsZeroed = false; - - for (const Cell &cell : column) { - auto it1 = column_.find(cell); - if (it1 != column_.end()) { - it1->get_element() = operators_->multiply_and_add(cell.get_element(), val, it1->get_element()); - if (it1->get_element() == Field_operators::get_additive_identity()){ - if constexpr (Master_matrix::isNonBasic && !Master_matrix::Option_list::is_of_boundary_type){ - if (it1->get_row_index() == chain_opt::get_pivot()) pivotIsZeroed = true; - } - _delete_cell(it1); - } else { - if constexpr (Master_matrix::Option_list::has_row_access) - ra_opt::update_cell(*it1); - } - } else { - _insert_cell(operators_->multiply(cell.get_element(), val), cell.get_row_index(), column_.end()); - } - } - - return pivotIsZeroed; -} - -} //namespace persistence_matrix -} //namespace Gudhi - -template -struct std::hash > -{ - size_t operator()(const Gudhi::persistence_matrix::Intrusive_set_column& column) const - { - std::size_t seed = 0; - for (const auto& cell : column){ - seed ^= std::hash()(cell.get_row_index() * static_cast(cell.get_element())) + 0x9e3779b9 + (seed << 6) + (seed >> 2); - } - return seed; - } + if (val == Field_operators::get_multiplicative_identity()) return *this; + + for (Cell& cell : column_) { + cell.get_element() = operators_->multiply(cell.get_element(), val); + if constexpr (Master_matrix::Option_list::has_row_access) ra_opt::update_cell(cell); + } + } + + return *this; +} + +template +template +inline Intrusive_set_column& +Intrusive_set_column::multiply_and_add(const Field_element_type& val, + const Cell_range& column) +{ + static_assert((!Master_matrix::isNonBasic || std::is_same_v), + "For boundary columns, the range has to be a column of same type to help ensure the validity of the " + "base element."); // could be removed, if we give the responsability to the user. + static_assert((!Master_matrix::isNonBasic || Master_matrix::Option_list::is_of_boundary_type), + "For chain columns, the given column cannot be constant."); + + if constexpr (Master_matrix::Option_list::is_z2) { + if (val) { + _add(column); + } else { + clear(); + _add(column); + } + } else { + _multiply_and_add(val, column); + } + + return *this; +} + +template +inline Intrusive_set_column& +Intrusive_set_column::multiply_and_add(const Field_element_type& val, + Intrusive_set_column& column) +{ + if constexpr (Master_matrix::isNonBasic && !Master_matrix::Option_list::is_of_boundary_type) { + // assumes that the addition never zeros out this column. + if constexpr (Master_matrix::Option_list::is_z2) { + if (val) { + if (_add(column)) { + chain_opt::swap_pivots(column); + dim_opt::swap_dimension(column); + } + } else { + throw std::invalid_argument("A chain column should not be multiplied by 0."); + } + } else { + if (_multiply_and_add(val, column)) { + chain_opt::swap_pivots(column); + dim_opt::swap_dimension(column); + } + } + } else { + if constexpr (Master_matrix::Option_list::is_z2) { + if (val) { + _add(column); + } else { + clear(); + _add(column); + } + } else { + _multiply_and_add(val, column); + } + } + + return *this; +} + +template +template +inline Intrusive_set_column& +Intrusive_set_column::multiply_and_add(const Cell_range& column, + const Field_element_type& val) +{ + static_assert((!Master_matrix::isNonBasic || std::is_same_v), + "For boundary columns, the range has to be a column of same type to help ensure the validity of the " + "base element."); // could be removed, if we give the responsability to the user. + static_assert((!Master_matrix::isNonBasic || Master_matrix::Option_list::is_of_boundary_type), + "For chain columns, the given column cannot be constant."); + + if constexpr (Master_matrix::Option_list::is_z2) { + if (val) { + _add(column); + } + } else { + _multiply_and_add(column, val); + } + + return *this; +} + +template +inline Intrusive_set_column& +Intrusive_set_column::multiply_and_add(Intrusive_set_column& column, + const Field_element_type& val) +{ + if constexpr (Master_matrix::isNonBasic && !Master_matrix::Option_list::is_of_boundary_type) { + // assumes that the addition never zeros out this column. + if constexpr (Master_matrix::Option_list::is_z2) { + if (val) { + if (_add(column)) { + chain_opt::swap_pivots(column); + dim_opt::swap_dimension(column); + } + } + } else { + if (_multiply_and_add(column, val)) { + chain_opt::swap_pivots(column); + dim_opt::swap_dimension(column); + } + } + } else { + if constexpr (Master_matrix::Option_list::is_z2) { + if (val) { + _add(column); + } + } else { + _multiply_and_add(column, val); + } + } + + return *this; +} + +template +inline Intrusive_set_column& +Intrusive_set_column::operator=(const Intrusive_set_column& other) +{ + static_assert(!Master_matrix::Option_list::has_row_access, "= assignement not enabled with row access option."); + + dim_opt::operator=(other); + chain_opt::operator=(other); + + // order is important + column_.clear_and_dispose(delete_disposer(this)); + operators_ = other.operators_; + cellPool_ = other.cellPool_; + column_.clone_from(other.column_, new_cloner(cellPool_), delete_disposer(this)); + + return *this; +} + +template +inline void Intrusive_set_column::_delete_cell(iterator& it) +{ + it = column_.erase_and_dispose(it, delete_disposer(this)); +} + +template +inline void Intrusive_set_column::_insert_cell(const Field_element_type& value, + id_index rowIndex, + const iterator& position) +{ + if constexpr (Master_matrix::Option_list::has_row_access) { + Cell* new_cell = cellPool_->construct(ra_opt::columnIndex_, rowIndex); + new_cell->set_element(value); + column_.insert(position, *new_cell); + ra_opt::insert_cell(rowIndex, new_cell); + } else { + Cell* new_cell = cellPool_->construct(rowIndex); + new_cell->set_element(value); + column_.insert(position, *new_cell); + } +} + +template +inline void Intrusive_set_column::_insert_cell(id_index rowIndex, + const iterator& position) +{ + if constexpr (Master_matrix::Option_list::has_row_access) { + Cell* new_cell = cellPool_->construct(ra_opt::columnIndex_, rowIndex); + column_.insert(position, *new_cell); + ra_opt::insert_cell(rowIndex, new_cell); + } else { + Cell* new_cell = cellPool_->construct(rowIndex); + column_.insert(position, *new_cell); + } +} + +template +template +inline bool Intrusive_set_column::_add(const Cell_range& column) +{ + bool pivotIsZeroed = false; + + for (const Cell& cell : column) { + auto it1 = column_.find(cell); + if (it1 != column_.end()) { + if constexpr (Master_matrix::Option_list::is_z2) { + if constexpr (Master_matrix::isNonBasic && !Master_matrix::Option_list::is_of_boundary_type) { + if (it1->get_row_index() == chain_opt::get_pivot()) pivotIsZeroed = true; + } + _delete_cell(it1); + } else { + it1->get_element() = operators_->add(it1->get_element(), cell.get_element()); + if (it1->get_element() == Field_operators::get_additive_identity()) { + if constexpr (Master_matrix::isNonBasic && !Master_matrix::Option_list::is_of_boundary_type) { + if (it1->get_row_index() == chain_opt::get_pivot()) pivotIsZeroed = true; + } + _delete_cell(it1); + } else { + if constexpr (Master_matrix::Option_list::has_row_access) ra_opt::update_cell(*it1); + } + } + } else { + if constexpr (Master_matrix::Option_list::is_z2) { + _insert_cell(cell.get_row_index(), column_.end()); + } else { + _insert_cell(cell.get_element(), cell.get_row_index(), column_.end()); + } + } + } + + return pivotIsZeroed; +} + +template +template +inline bool Intrusive_set_column::_multiply_and_add(const Field_element_type& val, + const Cell_range& column) +{ + bool pivotIsZeroed = false; + + if (val == 0u) { + if constexpr (Master_matrix::isNonBasic && !Master_matrix::Option_list::is_of_boundary_type) { + throw std::invalid_argument("A chain column should not be multiplied by 0."); + // this would not only mess up the base, but also the pivots stored. + } else { + clear(); + } + } + + // cannot use usual addition method, because all values of column_ have to be multiplied + + auto itTarget = column_.begin(); + auto itSource = column.begin(); + while (itTarget != column_.end() && itSource != column.end()) { + if (itTarget->get_row_index() < itSource->get_row_index()) { + itTarget->get_element() = operators_->multiply(itTarget->get_element(), val); + if constexpr (Master_matrix::Option_list::has_row_access) ra_opt::update_cell(*itTarget); + ++itTarget; + } else if (itTarget->get_row_index() > itSource->get_row_index()) { + _insert_cell(itSource->get_element(), itSource->get_row_index(), itTarget); + ++itSource; + } else { + itTarget->get_element() = operators_->multiply_and_add(itTarget->get_element(), val, itSource->get_element()); + if (itTarget->get_element() == Field_operators::get_additive_identity()) { + if constexpr (Master_matrix::isNonBasic && !Master_matrix::Option_list::is_of_boundary_type) { + if (itTarget->get_row_index() == chain_opt::get_pivot()) pivotIsZeroed = true; + } + _delete_cell(itTarget); + } else { + if constexpr (Master_matrix::Option_list::has_row_access) ra_opt::update_cell(*itTarget); + ++itTarget; + } + ++itSource; + } + } + + while (itTarget != column_.end()) { + itTarget->get_element() = operators_->multiply(itTarget->get_element(), val); + if constexpr (Master_matrix::Option_list::has_row_access) ra_opt::update_cell(*itTarget); + itTarget++; + } + + while (itSource != column.end()) { + _insert_cell(itSource->get_element(), itSource->get_row_index(), column_.end()); + ++itSource; + } + + return pivotIsZeroed; +} + +template +template +inline bool Intrusive_set_column::_multiply_and_add(const Cell_range& column, + const Field_element_type& val) +{ + if (val == 0u) { + return false; + } + + bool pivotIsZeroed = false; + + for (const Cell& cell : column) { + auto it1 = column_.find(cell); + if (it1 != column_.end()) { + it1->get_element() = operators_->multiply_and_add(cell.get_element(), val, it1->get_element()); + if (it1->get_element() == Field_operators::get_additive_identity()) { + if constexpr (Master_matrix::isNonBasic && !Master_matrix::Option_list::is_of_boundary_type) { + if (it1->get_row_index() == chain_opt::get_pivot()) pivotIsZeroed = true; + } + _delete_cell(it1); + } else { + if constexpr (Master_matrix::Option_list::has_row_access) ra_opt::update_cell(*it1); + } + } else { + _insert_cell(operators_->multiply(cell.get_element(), val), cell.get_row_index(), column_.end()); + } + } + + return pivotIsZeroed; +} + +} // namespace persistence_matrix +} // namespace Gudhi + + +/** + * @brief Hash method for @ref Gudhi::persistence_matrix::Intrusive_set_column. + * + * @tparam Master_matrix Template parameter of @ref Gudhi::persistence_matrix::Intrusive_set_column. + * @tparam Cell_constructor Template parameter of @ref Gudhi::persistence_matrix::Intrusive_set_column. + */ +template +struct std::hash > +{ + size_t operator()( + const Gudhi::persistence_matrix::Intrusive_set_column& column) const { + std::size_t seed = 0; + for (const auto& cell : column) { + seed ^= std::hash()(cell.get_row_index() * static_cast(cell.get_element())) + + 0x9e3779b9 + (seed << 6) + (seed >> 2); + } + return seed; + } }; -#endif // PM_INTRUSIVE_SET_COLUMN_H +#endif // PM_INTRUSIVE_SET_COLUMN_H diff --git a/src/Persistence_matrix/include/gudhi/matrix.h b/src/Persistence_matrix/include/gudhi/matrix.h index feedf3097d..c704cb0339 100644 --- a/src/Persistence_matrix/include/gudhi/matrix.h +++ b/src/Persistence_matrix/include/gudhi/matrix.h @@ -126,7 +126,7 @@ class Matrix { using index = typename PersistenceMatrixOptions::index_type; /**< Type of MatIdx index. */ using id_index = typename PersistenceMatrixOptions::index_type; /**< Type of IDIdx index. */ using pos_index = typename PersistenceMatrixOptions::index_type; /**< Type of PosIdx index. */ - using dimension_type = typename PersistenceMatrixOptions::dimension_type; /**< Type for dimension. */ + using dimension_type = typename PersistenceMatrixOptions::dimension_type; /**< Type for dimension value. */ struct Dummy_field_operators{ using element_type = unsigned int; @@ -338,7 +338,8 @@ class Matrix { /** * @brief Type of the columns stored in the matrix. The type depends on the value of @ref column_type defined - * in the given options. See @ref Column_types for a more detailed description. + * in the given options. See @ref Column_types for a more detailed description. All columns follow the + * @ref PersistenceMatrixColumn concept. */ using Column_type = typename std::conditional< PersistenceMatrixOptions::column_type == Column_types::HEAP,