Matrix.h

Go to the documentation of this file.

/*    This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT.
 *    See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details.
 *    Author(s):       Hannah Schreiber
 *
 *    Copyright (C) 2022-24 Inria
 *
 *    Modification(s):
 *      - YYYY/MM Author: Description of the modification
 */
 
#ifndef MASTER_MATRIX_H
#define MASTER_MATRIX_H
 
#include <type_traits>
#include <vector>
#include <unordered_map>
#include <map>
#include <initializer_list>
 
#include <boost/intrusive/list.hpp>
 
#include <gudhi/Debug_utils.h>
 
#include <gudhi/persistence_matrix_options.h>
#include <gudhi/persistence_interval.h>
 
#include <gudhi/Fields/Z2_field_operators.h>
 
#include <gudhi/Persistence_matrix/Id_to_index_overlay.h>
#include <gudhi/Persistence_matrix/Position_to_index_overlay.h>
 
#include <gudhi/Persistence_matrix/matrix_dimension_holders.h>
#include <gudhi/Persistence_matrix/matrix_row_access.h>
#include <gudhi/Persistence_matrix/base_swap.h>
#include <gudhi/Persistence_matrix/base_pairing.h>
#include <gudhi/Persistence_matrix/ru_pairing.h>
#include <gudhi/Persistence_matrix/ru_vine_swap.h>
#include <gudhi/Persistence_matrix/ru_rep_cycles.h>
#include <gudhi/Persistence_matrix/chain_pairing.h>
#include <gudhi/Persistence_matrix/chain_vine_swap.h>
#include <gudhi/Persistence_matrix/chain_rep_cycles.h>
 
#include <gudhi/Persistence_matrix/Base_matrix.h>
#include <gudhi/Persistence_matrix/Base_matrix_with_column_compression.h>
#include <gudhi/Persistence_matrix/Boundary_matrix.h>
#include <gudhi/Persistence_matrix/RU_matrix.h>
#include <gudhi/Persistence_matrix/Chain_matrix.h>
 
#include <gudhi/Persistence_matrix/allocators/entry_constructors.h>
#include <gudhi/Persistence_matrix/columns/entry_types.h>
#include <gudhi/Persistence_matrix/columns/row_access.h>
 
#include <gudhi/Persistence_matrix/columns/column_dimension_holder.h>
#include <gudhi/Persistence_matrix/columns/chain_column_extra_properties.h>
#include <gudhi/Persistence_matrix/columns/intrusive_list_column.h>
#include <gudhi/Persistence_matrix/columns/intrusive_set_column.h>
#include <gudhi/Persistence_matrix/columns/list_column.h>
#include <gudhi/Persistence_matrix/columns/set_column.h>
#include <gudhi/Persistence_matrix/columns/unordered_set_column.h>
#include <gudhi/Persistence_matrix/columns/vector_column.h>
#include <gudhi/Persistence_matrix/columns/naive_vector_column.h>
#include <gudhi/Persistence_matrix/columns/heap_column.h>
 
namespace Gudhi {
namespace persistence_matrix {
 
template <class PersistenceMatrixOptions = Default_options<> >
class Matrix {
 public:
  using Option_list = PersistenceMatrixOptions;     //to make it accessible from the other classes
  using Index = typename PersistenceMatrixOptions::Index;         
  using ID_index = typename PersistenceMatrixOptions::Index;      
  using Pos_index = typename PersistenceMatrixOptions::Index;     
  using Dimension = typename PersistenceMatrixOptions::Dimension; 
  using Field_operators =
      typename std::conditional<PersistenceMatrixOptions::is_z2,
                                Gudhi::persistence_fields::Z2_field_operators,
                                typename PersistenceMatrixOptions::Field_coeff_operators
                               >::type;
  using Element = typename Field_operators::Element;
  using Characteristic = typename Field_operators::Characteristic;
 
  template <typename T>
  static constexpr const T get_null_value() {
    return -1;
  }
 
  using Bar = Persistence_interval<Dimension, Pos_index>;
 
  //tags for boost to associate a row and a column to a same entry
  struct Matrix_row_tag;
  struct Matrix_column_tag;
 
  using Base_hook_matrix_row =
      boost::intrusive::list_base_hook<boost::intrusive::tag<Matrix_row_tag>,
                                       boost::intrusive::link_mode<boost::intrusive::auto_unlink> >;
  using Base_hook_matrix_list_column =
      boost::intrusive::list_base_hook<boost::intrusive::tag<Matrix_column_tag>,
                                       boost::intrusive::link_mode<boost::intrusive::safe_link> >;
  using Base_hook_matrix_set_column =
      boost::intrusive::set_base_hook<boost::intrusive::tag<Matrix_column_tag>,
                                      boost::intrusive::link_mode<boost::intrusive::safe_link> >;
 
  //Two dummies are necessary to avoid double inheritance as an entry can inherit both a row and a column hook.
  struct Dummy_row_hook {};
  struct Dummy_column_hook {};
 
  using Row_hook = typename std::conditional<PersistenceMatrixOptions::has_row_access &&
                                                 PersistenceMatrixOptions::has_intrusive_rows,
                                             Base_hook_matrix_row,
                                             Dummy_row_hook
                                            >::type;
  using Column_hook = typename std::conditional<
      PersistenceMatrixOptions::column_type == Column_types::INTRUSIVE_LIST,
      Base_hook_matrix_list_column,
      typename std::conditional<PersistenceMatrixOptions::column_type == Column_types::INTRUSIVE_SET,
                                Base_hook_matrix_set_column,
                                Dummy_column_hook
                               >::type
    >::type;
 
  //Option to store the column index within the entry (additionally to the row index). Necessary only with row access.
  using Entry_column_index_option =
      typename std::conditional<PersistenceMatrixOptions::has_row_access,
                                Entry_column_index<Index>,
                                Dummy_entry_column_index_mixin
                               >::type;
  //Option to store the value of the entry.
  //Unnecessary for values in Z_2 as there are always 1 (0-valued entries are never stored).
  using Entry_field_element_option =
      typename std::conditional<PersistenceMatrixOptions::is_z2,
                                Dummy_entry_field_element_mixin,
                                Entry_field_element<Element>
                               >::type;
  using Matrix_entry = Entry<Matrix<PersistenceMatrixOptions> >;
 
  inline static New_entry_constructor<Matrix_entry> defaultEntryConstructor;
  using Entry_constructor = Pool_entry_constructor<Matrix_entry>;
 
  using Entry_representative = typename std::conditional<PersistenceMatrixOptions::is_z2,
                                                  ID_index,
                                                  std::pair<ID_index, Element>
                                                 >::type;
 
  struct RowEntryComp {
    bool operator()(const Matrix_entry& c1, const Matrix_entry& c2) const {
      return c1.get_column_index() < c2.get_column_index();
    }
  };
 
  using Row =
      typename std::conditional<PersistenceMatrixOptions::has_intrusive_rows,
                                boost::intrusive::list<Matrix_entry,
                                                       boost::intrusive::constant_time_size<false>,
                                                       boost::intrusive::base_hook<Base_hook_matrix_row>
                                                      >,
                                std::set<Matrix_entry, RowEntryComp>
                               >::type;
 
  using Row_container =
      typename std::conditional<PersistenceMatrixOptions::has_removable_rows,
                                std::map<ID_index, Row>,
                                std::vector<Row>
                               >::type;
 
  //Row access at column level
  using Row_access_option =
      typename std::conditional<PersistenceMatrixOptions::has_row_access,
                                Row_access<Matrix<PersistenceMatrixOptions> >,
                                Dummy_row_access
                               >::type;
  //Row access at matrix level
  using Matrix_row_access_option =
      typename std::conditional<PersistenceMatrixOptions::has_row_access,
                                Matrix_row_access<Row,
                                                  Row_container,
                                                  PersistenceMatrixOptions::has_removable_rows,
                                                  ID_index>,
                                Dummy_matrix_row_access
                               >::type;
 
  template <typename value_type>
  using Dictionary =
      typename std::conditional<PersistenceMatrixOptions::has_map_column_container,
                                std::unordered_map<unsigned int, value_type>,
                                std::vector<value_type>
                               >::type;
 
  static const bool isNonBasic = PersistenceMatrixOptions::has_column_pairings ||
                                 PersistenceMatrixOptions::has_vine_update ||
                                 PersistenceMatrixOptions::can_retrieve_representative_cycles;
 
  using Column_dimension_option =
      typename std::conditional<isNonBasic,
                                Column_dimension_holder<Matrix<PersistenceMatrixOptions> >,
                                Dummy_dimension_holder
                               >::type;
  //Extra information needed for a column when the matrix is a @ref chainmatrix "chain matrix".
  using Chain_column_option =
      typename std::conditional<isNonBasic && !PersistenceMatrixOptions::is_of_boundary_type,
                                Chain_column_extra_properties<Matrix<PersistenceMatrixOptions> >,
                                Dummy_chain_properties
                               >::type;
 
  using Matrix_heap_column = Heap_column<Matrix<PersistenceMatrixOptions> >;
  using Matrix_list_column = List_column<Matrix<PersistenceMatrixOptions> >;
  using Matrix_vector_column = Vector_column<Matrix<PersistenceMatrixOptions> >;
  using Matrix_naive_vector_column = Naive_std_vector_column<Matrix<PersistenceMatrixOptions> >;
  using Matrix_small_vector_column = Naive_small_vector_column<Matrix<PersistenceMatrixOptions> >;
  using Matrix_set_column = Set_column<Matrix<PersistenceMatrixOptions> >;
  using Matrix_unordered_set_column = Unordered_set_column<Matrix<PersistenceMatrixOptions> >;
  using Matrix_intrusive_list_column = Intrusive_list_column<Matrix<PersistenceMatrixOptions> >;
  using Matrix_intrusive_set_column = Intrusive_set_column<Matrix<PersistenceMatrixOptions> >;
 
  using Column = typename std::conditional<
        PersistenceMatrixOptions::column_type == Column_types::HEAP,
        Matrix_heap_column,
        typename std::conditional<
            PersistenceMatrixOptions::column_type == Column_types::LIST,
            Matrix_list_column,
            typename std::conditional<
                PersistenceMatrixOptions::column_type == Column_types::SET,
                Matrix_set_column,
                typename std::conditional<
                    PersistenceMatrixOptions::column_type == Column_types::UNORDERED_SET,
                    Matrix_unordered_set_column,
                    typename std::conditional<
                        PersistenceMatrixOptions::column_type == Column_types::VECTOR,
                        Matrix_vector_column,
                        typename std::conditional<
                            PersistenceMatrixOptions::column_type == Column_types::INTRUSIVE_LIST,
                            Matrix_intrusive_list_column,
                            typename std::conditional<
                                PersistenceMatrixOptions::column_type == Column_types::NAIVE_VECTOR,
                                Matrix_naive_vector_column,
                                typename std::conditional<
                                    PersistenceMatrixOptions::column_type == Column_types::SMALL_VECTOR,
                                    Matrix_small_vector_column,
                                    Matrix_intrusive_set_column
                                    >::type
                                >::type
                            >::type
                        >::type
                    >::type
                >::type
            >::type
        >::type;
 
  struct Column_z2_settings{
    Column_z2_settings() : entryConstructor() {}
    Column_z2_settings([[maybe_unused]] Characteristic characteristic) : entryConstructor() {}
    Column_z2_settings(const Column_z2_settings& toCopy) : entryConstructor() {}
 
    Entry_constructor entryConstructor;   //will be replaced by more specific allocators depending on the column type.
  };
 
  struct Column_zp_settings {
    Column_zp_settings() : operators(), entryConstructor() {}
    //purposely triggers operators() instead of operators(characteristic) as the "dummy" values for the different
    //operators can be different from -1.
    Column_zp_settings(Characteristic characteristic) : operators(), entryConstructor() {
      if (characteristic != get_null_value<Characteristic>()) operators.set_characteristic(characteristic);
    }
    Column_zp_settings(const Column_zp_settings& toCopy)
        : operators(toCopy.operators.get_characteristic()), entryConstructor() {}
 
    Field_operators operators;
    Entry_constructor entryConstructor;   //will be replaced by more specific allocators depending on the column type.
  };
 
  // struct Column_z2_with_rows_settings {
  //   Column_z2_with_rows_settings() : entryConstructor(), rows(nullptr) {}
  //   Column_z2_with_rows_settings([[maybe_unused]] Characteristic characteristic)
  //       : entryConstructor(), rows(nullptr) {}
 
  //   Entry_constructor entryConstructor;
  //   Row_container* rows;
  // };
 
  // struct Column_zp_with_rows_settings {
  //   Column_zp_with_rows_settings() : operators(), entryConstructor(), rows(nullptr) {}
  //   Column_zp_with_rows_settings(Characteristic characteristic)
  //       : operators(characteristic), entryConstructor(), rows(nullptr) {}
 
  //   Field_operators operators;
  //   Entry_constructor entryConstructor;
  //   Row_container* rows;
  // };
 
  //To prepare a more flexible use of the column types later (custom allocators depending on the column type etc.)
  using Column_settings = typename std::conditional<
      PersistenceMatrixOptions::is_z2,
      Column_z2_settings,
      Column_zp_settings
    >::type;
 
  // using Column_settings = typename std::conditional<
  //     PersistenceMatrixOptions::is_z2,
  //     typename std::conditional<PersistenceMatrixOptions::has_row_access,
  //                               Column_z2_with_rows_settings,
  //                               Column_z2_settings
  //                              >::type,
  //     typename std::conditional<PersistenceMatrixOptions::has_row_access,
  //                               Column_zp_with_rows_settings,
  //                               Column_zp_settings
  //                              >::type
  //   >::type;
 
  using Column_container =
      typename std::conditional<PersistenceMatrixOptions::has_map_column_container,
                                std::unordered_map<Index, Column>,
                                std::vector<Column>
                               >::type;
 
  static const bool hasFixedBarcode = Option_list::is_of_boundary_type && !PersistenceMatrixOptions::has_vine_update;
  using Barcode = typename std::conditional<hasFixedBarcode,
                                            std::vector<Bar>,
                                            typename std::conditional<PersistenceMatrixOptions::has_removable_columns,
                                                                      std::list<Bar>,
                                                                      std::vector<Bar>
                                                                     >::type
                                           >::type;
  using Bar_dictionary =
      typename std::conditional<hasFixedBarcode,
                                typename std::conditional<PersistenceMatrixOptions::can_retrieve_representative_cycles,
                                  std::vector<Index>,                   //RU
                                  std::unordered_map<Pos_index, Index>  //boundary
                                >::type,
                                typename std::conditional<PersistenceMatrixOptions::has_removable_columns,
                                  std::unordered_map<Pos_index, typename Barcode::iterator>,
                                  std::vector<Index>
                                >::type
                               >::type;
 
  //default type for boundaries to permit list initialization directly in function parameters
  using Boundary = typename std::conditional<PersistenceMatrixOptions::is_z2,
                                             std::initializer_list<ID_index>,
                                             std::initializer_list<std::pair<ID_index, Element> >
                                            >::type;
 
  //i.e. is simple @ref boundarymatrix "boundary matrix". Also, only needed because of the reduction algorithm.
  //TODO: remove the necessity and recalculate when needed or keep it like that?
  static const bool maxDimensionIsNeeded =
      PersistenceMatrixOptions::has_column_pairings && PersistenceMatrixOptions::is_of_boundary_type &&
      !PersistenceMatrixOptions::has_vine_update && !PersistenceMatrixOptions::can_retrieve_representative_cycles;
 
  using Matrix_dimension_option = typename std::conditional<
      PersistenceMatrixOptions::has_matrix_maximal_dimension_access || maxDimensionIsNeeded,
      typename std::conditional<PersistenceMatrixOptions::has_removable_columns,
                                Matrix_all_dimension_holder<Dimension>,
                                Matrix_max_dimension_holder<Dimension>
                               >::type,
      Dummy_matrix_dimension_holder
    >::type;
 
  using Master_base_matrix =
      typename std::conditional<PersistenceMatrixOptions::has_column_compression,
                                Base_matrix_with_column_compression<Matrix<PersistenceMatrixOptions> >,
                                Base_matrix<Matrix<PersistenceMatrixOptions> >
                               >::type;
  using Master_boundary_matrix = Boundary_matrix<Matrix<PersistenceMatrixOptions> >;
  using Master_RU_matrix = RU_matrix<Matrix<PersistenceMatrixOptions> >;
  using Master_chain_matrix = Chain_matrix<Matrix<PersistenceMatrixOptions> >;
 
  template <class Base>
  using Base_swap_option =
      typename std::conditional<PersistenceMatrixOptions::has_vine_update ||
                                    PersistenceMatrixOptions::has_column_and_row_swaps,
                                Base_swap<Matrix<PersistenceMatrixOptions>, Base>,
                                Dummy_base_swap
                               >::type;
  using Base_pairing_option =
      typename std::conditional<PersistenceMatrixOptions::has_column_pairings &&
                                    !PersistenceMatrixOptions::has_vine_update &&
                                    !PersistenceMatrixOptions::can_retrieve_representative_cycles,
                                Base_pairing<Matrix<PersistenceMatrixOptions> >,
                                Dummy_base_pairing
                               >::type;
 
  using RU_pairing_option =
      typename std::conditional<PersistenceMatrixOptions::has_column_pairings &&
                                    !PersistenceMatrixOptions::has_vine_update,
                                RU_pairing<Matrix<PersistenceMatrixOptions> >,
                                Dummy_ru_pairing
                               >::type;
  using RU_vine_swap_option =
      typename std::conditional<PersistenceMatrixOptions::has_vine_update,
                                RU_vine_swap<Matrix<PersistenceMatrixOptions> >,
                                Dummy_ru_vine_swap
                               >::type;
  using RU_representative_cycles_option =
      typename std::conditional<PersistenceMatrixOptions::can_retrieve_representative_cycles,
                                RU_representative_cycles<Matrix<PersistenceMatrixOptions> >,
                                Dummy_ru_representative_cycles
                               >::type;
 
  using Chain_pairing_option =
      typename std::conditional<PersistenceMatrixOptions::has_column_pairings &&
                                    !PersistenceMatrixOptions::has_vine_update,
                                Chain_pairing<Matrix<PersistenceMatrixOptions> >,
                                Dummy_chain_pairing
                               >::type;
  using Chain_vine_swap_option = typename std::conditional<PersistenceMatrixOptions::has_vine_update,
                                                           Chain_vine_swap<Matrix<PersistenceMatrixOptions> >,
                                                           Dummy_chain_vine_swap
                                                          >::type;
  using Chain_representative_cycles_option =
      typename std::conditional<PersistenceMatrixOptions::can_retrieve_representative_cycles,
                                Chain_representative_cycles<Matrix<PersistenceMatrixOptions> >,
                                Dummy_chain_representative_cycles
                               >::type;
 
  using Cycle = std::vector<ID_index>; //TODO: add coefficients
 
  //Return types to factorize the corresponding methods
 
  //The returned column is `const` if the matrix uses column compression
  using Returned_column =
      typename std::conditional<!isNonBasic && PersistenceMatrixOptions::has_column_compression,
                                const Column,
                                Column
                               >::type;
  //The returned row is `const` if the matrix uses column compression
  using Returned_row =
      typename std::conditional<!isNonBasic && PersistenceMatrixOptions::has_column_compression,
                                const Row,
                                Row
                               >::type;
  //If the matrix is a chain matrix, the insertion method returns the pivots of its unpaired columns used to reduce the
  //inserted boundary. Otherwise, void.
  using Insertion_return =
      typename std::conditional<PersistenceMatrixOptions::is_of_boundary_type || !isNonBasic ||
                                    PersistenceMatrixOptions::column_indexation_type ==
                                        Column_indexation_types::POSITION,
                                void,
                                std::vector<Entry_representative>
                               >::type;
 
  Matrix();
  template <class Container = Boundary>
  Matrix(const std::vector<Container>& columns, Characteristic characteristic = 11);
  Matrix(unsigned int numberOfColumns, Characteristic characteristic = Matrix::get_null_value<Characteristic>());
  Matrix(const std::function<bool(Pos_index,Pos_index)>& birthComparator,
         const std::function<bool(Pos_index,Pos_index)>& deathComparator);
  template <class Boundary_range = Boundary>
  Matrix(const std::vector<Boundary_range>& orderedBoundaries,
         const std::function<bool(Pos_index,Pos_index)>& birthComparator,
         const std::function<bool(Pos_index,Pos_index)>& deathComparator,
         Characteristic characteristic = 11);
  Matrix(unsigned int numberOfColumns,
         const std::function<bool(Pos_index,Pos_index)>& birthComparator,
         const std::function<bool(Pos_index,Pos_index)>& deathComparator,
         Characteristic characteristic = Matrix::get_null_value<Characteristic>());
  Matrix(const Matrix& matrixToCopy);
  Matrix(Matrix&& other) noexcept;
 
  ~Matrix();
 
  //TODO: compatibility with multi fields:
  //  - set_characteristic(Characteristic min, Characteristic max)
  //  - readapt reduction?
  void set_characteristic(Characteristic characteristic);
 
  // (TODO: if there is no row access and the column type corresponds to the internal column type of the matrix,
  // moving the column instead of copying it should be possible. Is it worth implementing it?)
  template <class Container>
  void insert_column(const Container& column);
  template <class Container>
  void insert_column(const Container& column, Index columnIndex);
  //TODO: for simple boundary matrices, add an index pointing to the first column inserted after the last call of
  //get_current_barcode to enable several calls to get_current_barcode
  template <class Boundary_range = Boundary>
  Insertion_return insert_boundary(const Boundary_range& boundary, Dimension dim = Matrix::get_null_value<Dimension>());
  template <class Boundary_range = Boundary>
  Insertion_return insert_boundary(ID_index cellIndex, const Boundary_range& boundary,
                                   Dimension dim = Matrix::get_null_value<Dimension>());
 
  Returned_column& get_column(Index columnIndex);
  const Column& get_column(Index columnIndex) const;
  //TODO: there is no particular reason that this method is not available for identifier indexing,
  // just has to be added to the interface...
  const Column& get_column(Index columnIndex, bool inR);
 
  //TODO: update column indices when reordering rows (after lazy swap) such that always MatIdx are returned.
  Returned_row& get_row(ID_index rowIndex);
  const Row& get_row(ID_index rowIndex) const;
  //TODO: there is no particular reason that this method is not available for identifier indexing,
  // just has to be added to the interface...
  const Row& get_row(ID_index rowIndex, bool inR);
 
  void remove_column(Index columnIndex);
  //TODO: rename method to be less confusing.
  void erase_empty_row(ID_index rowIndex);
  //TODO: for chain matrices, replace IDIdx input with MatIdx input to homogenize.
  void remove_maximal_cell(Index columnIndex);
  //TODO: See if it would be better to use something more general than a vector for columnsToSwap, such that
  // the user do not have to construct the vector from scratch. Like passing iterators instead. But it would be nice,
  // to still be able to do (cell, {})...
  void remove_maximal_cell(ID_index cellIndex, const std::vector<ID_index>& columnsToSwap);
  void remove_last();
 
  Dimension get_max_dimension() const;
  Index get_number_of_columns() const;
  Dimension get_column_dimension(Index columnIndex) const;
 
  template <typename Integer_index>
  std::enable_if_t<std::is_integral_v<Integer_index> > add_to(Integer_index sourceColumnIndex,
                                                              Integer_index targetColumnIndex);
  template <class Entry_range>
  std::enable_if_t<!std::is_integral_v<Entry_range> > add_to(const Entry_range& sourceColumn, Index targetColumnIndex);
 
  template <typename Integer_index>
  std::enable_if_t<std::is_integral_v<Integer_index> > multiply_target_and_add_to(Integer_index sourceColumnIndex,
                                                                               int coefficient,
                                                                               Integer_index targetColumnIndex);
  template <class Entry_range>
  std::enable_if_t<!std::is_integral_v<Entry_range> > multiply_target_and_add_to(const Entry_range& sourceColumn,
                                                                                int coefficient,
                                                                                Index targetColumnIndex);
 
  template <typename Integer_index>
  std::enable_if_t<std::is_integral_v<Integer_index> > multiply_source_and_add_to(int coefficient,
                                                                               Integer_index sourceColumnIndex,
                                                                               Integer_index targetColumnIndex);
  template <class Entry_range>
  std::enable_if_t<!std::is_integral_v<Entry_range> > multiply_source_and_add_to(int coefficient,
                                                                                const Entry_range& sourceColumn,
                                                                                Index targetColumnIndex);
 
  void zero_entry(Index columnIndex, ID_index rowIndex);
  void zero_entry(Index columnIndex, ID_index rowIndex, bool inR);
  void zero_column(Index columnIndex);
  void zero_column(Index columnIndex, bool inR);
  bool is_zero_entry(Index columnIndex, ID_index rowIndex);
  bool is_zero_entry(Index columnIndex, ID_index rowIndex, bool inR) const;
  bool is_zero_column(Index columnIndex);
  bool is_zero_column(Index columnIndex, bool inR);
 
  Index get_column_with_pivot(ID_index cellIndex) const;
  ID_index get_pivot(Index columnIndex);
 
  Matrix& operator=(Matrix other);
  friend void swap(Matrix& matrix1, Matrix& matrix2) {
    swap(matrix1.matrix_, matrix2.matrix_);
    std::swap(matrix1.colSettings_, matrix2.colSettings_);
  }
 
  void print();  // for debug
 
  //TODO: change the behaviour for boundary matrices.
  const Barcode& get_current_barcode();
  const Barcode& get_current_barcode() const;
 
  void swap_columns(Index columnIndex1, Index columnIndex2);
  void swap_rows(Index rowIndex1, Index rowIndex2);
  //TODO: find better name. And benchmark also to verify if it is really worth it to have this extra version in addition
  //to vine_swap.
  bool vine_swap_with_z_eq_1_case(Pos_index index);
  Index vine_swap_with_z_eq_1_case(Index columnIndex1, Index columnIndex2);
  bool vine_swap(Pos_index index);
  Index vine_swap(Index columnIndex1, Index columnIndex2);
 
  //TODO: Rethink the interface for representative cycles
  void update_representative_cycles();
  const std::vector<Cycle>& get_representative_cycles();
  const Cycle& get_representative_cycle(const Bar& bar);
 
 private:
  using Underlying_matrix =
    typename std::conditional<
        isNonBasic,
        typename std::conditional<
            PersistenceMatrixOptions::is_of_boundary_type,
            typename std::conditional<
                PersistenceMatrixOptions::has_vine_update ||
                  PersistenceMatrixOptions::can_retrieve_representative_cycles,
                typename std::conditional<
                    PersistenceMatrixOptions::column_indexation_type == Column_indexation_types::CONTAINER ||
                      PersistenceMatrixOptions::column_indexation_type == Column_indexation_types::POSITION,
                    Master_RU_matrix, Id_to_index_overlay<Master_RU_matrix, Matrix<PersistenceMatrixOptions> >
                >::type,
                typename std::conditional<
                    PersistenceMatrixOptions::column_indexation_type == Column_indexation_types::CONTAINER ||
                      PersistenceMatrixOptions::column_indexation_type == Column_indexation_types::POSITION,
                    Master_boundary_matrix,
                    Id_to_index_overlay<Master_boundary_matrix, Matrix<PersistenceMatrixOptions> >
                >::type
            >::type,
            typename std::conditional<
                PersistenceMatrixOptions::column_indexation_type == Column_indexation_types::CONTAINER,
                Master_chain_matrix,
                typename std::conditional<
                    PersistenceMatrixOptions::column_indexation_type == Column_indexation_types::POSITION,
                    Position_to_index_overlay<Master_chain_matrix, Matrix<PersistenceMatrixOptions> >,
                    Id_to_index_overlay<Master_chain_matrix, Matrix<PersistenceMatrixOptions> >
                >::type
            >::type
        >::type,
        Master_base_matrix
    >::type;
 
  // Field_operators* operators_;
  // Entry_constructor* entryPool_;
  Column_settings* colSettings_;  //pointer because the of swap operator on matrix_ which also stores the pointer
  Underlying_matrix matrix_;
 
  static constexpr void _assert_options();
};
 
template <class PersistenceMatrixOptions>
inline Matrix<PersistenceMatrixOptions>::Matrix()
    : colSettings_(new Column_settings()), matrix_(colSettings_)
{
  static_assert(
      PersistenceMatrixOptions::is_of_boundary_type || !PersistenceMatrixOptions::has_vine_update ||
          PersistenceMatrixOptions::has_column_pairings,
      "When no barcode is recorded with vine swaps, comparison functions for the columns have to be provided.");
  _assert_options();
}
 
template <class PersistenceMatrixOptions>
template <class Container>
inline Matrix<PersistenceMatrixOptions>::Matrix(const std::vector<Container>& columns,
                                                Characteristic characteristic)
    : colSettings_(new Column_settings(characteristic)),
      matrix_(columns, colSettings_)
{
  static_assert(PersistenceMatrixOptions::is_of_boundary_type || !PersistenceMatrixOptions::has_vine_update ||
                    PersistenceMatrixOptions::has_column_pairings,
                "When no barcode is recorded with vine swaps for chain matrices, comparison functions for the columns "
                "have to be provided.");
  _assert_options();
}
 
template <class PersistenceMatrixOptions>
inline Matrix<PersistenceMatrixOptions>::Matrix(unsigned int numberOfColumns, Characteristic characteristic)
    : colSettings_(new Column_settings(characteristic)),
      matrix_(numberOfColumns, colSettings_)
{
  static_assert(PersistenceMatrixOptions::is_of_boundary_type || !PersistenceMatrixOptions::has_vine_update ||
                    PersistenceMatrixOptions::has_column_pairings,
                "When no barcode is recorded with vine swaps for chain matrices, comparison functions for the columns "
                "have to be provided.");
  _assert_options();
}
 
template <class PersistenceMatrixOptions>
inline Matrix<PersistenceMatrixOptions>::Matrix(const std::function<bool(Pos_index,Pos_index)>& birthComparator,
                                                const std::function<bool(Pos_index,Pos_index)>& deathComparator)
    : colSettings_(new Column_settings()),
      matrix_(colSettings_, birthComparator, deathComparator)
{
  static_assert(
      !PersistenceMatrixOptions::is_of_boundary_type && PersistenceMatrixOptions::has_vine_update &&
          !PersistenceMatrixOptions::has_column_pairings,
      "Constructor only available for chain matrices when vine swaps are enabled, but the barcode is not recorded.");
  _assert_options();
}
 
template <class PersistenceMatrixOptions>
template <class Boundary_range>
inline Matrix<PersistenceMatrixOptions>::Matrix(const std::vector<Boundary_range>& orderedBoundaries,
                                                const std::function<bool(Pos_index, Pos_index)>& birthComparator,
                                                const std::function<bool(Pos_index, Pos_index)>& deathComparator,
                                                Characteristic characteristic)
    : colSettings_(new Column_settings(characteristic)),
      matrix_(orderedBoundaries, colSettings_, birthComparator, deathComparator)
{
  static_assert(
      !PersistenceMatrixOptions::is_of_boundary_type && PersistenceMatrixOptions::has_vine_update &&
          !PersistenceMatrixOptions::has_column_pairings,
      "Constructor only available for chain matrices when vine swaps are enabled, but the barcode is not recorded.");
  _assert_options();
}
 
template <class PersistenceMatrixOptions>
inline Matrix<PersistenceMatrixOptions>::Matrix(unsigned int numberOfColumns,
                                                const std::function<bool(Pos_index, Pos_index)>& birthComparator,
                                                const std::function<bool(Pos_index, Pos_index)>& deathComparator,
                                                Characteristic characteristic)
    : colSettings_(new Column_settings(characteristic)),
      matrix_(numberOfColumns, colSettings_, birthComparator, deathComparator)
{
  static_assert(
      !PersistenceMatrixOptions::is_of_boundary_type && PersistenceMatrixOptions::has_vine_update &&
          !PersistenceMatrixOptions::has_column_pairings,
      "Constructor only available for chain matrices when vine swaps are enabled, but the barcode is not recorded.");
  _assert_options();
}
 
template <class PersistenceMatrixOptions>
inline Matrix<PersistenceMatrixOptions>::Matrix(const Matrix& matrixToCopy)
    : colSettings_(new Column_settings(*matrixToCopy.colSettings_)),
      matrix_(matrixToCopy.matrix_, colSettings_)
{
  _assert_options();
}
 
template <class PersistenceMatrixOptions>
inline Matrix<PersistenceMatrixOptions>::Matrix(Matrix&& other) noexcept
    : colSettings_(std::exchange(other.colSettings_, nullptr)),
      matrix_(std::move(other.matrix_))
{
  _assert_options();
}
 
template <class PersistenceMatrixOptions>
inline Matrix<PersistenceMatrixOptions>::~Matrix()
{
  matrix_.reset(colSettings_);
  delete colSettings_;
}
 
template <class PersistenceMatrixOptions>
inline void Matrix<PersistenceMatrixOptions>::set_characteristic(Characteristic characteristic)
{
  if constexpr (!PersistenceMatrixOptions::is_z2) {
    if (colSettings_->operators.get_characteristic() != get_null_value<Characteristic>()) {
      std::cerr << "Warning: Characteristic already initialised. Changing it could lead to incoherences in the matrix "
                   "as the modulo was already applied to values in existing columns.";
    }
 
    colSettings_->operators.set_characteristic(characteristic);
  }
}
 
template <class PersistenceMatrixOptions>
template <class Container>
inline void Matrix<PersistenceMatrixOptions>::insert_column(const Container& column)
{
  if constexpr (!PersistenceMatrixOptions::is_z2){
    GUDHI_CHECK(colSettings_->operators.get_characteristic() != get_null_value<Characteristic>(),
                std::logic_error("Matrix::insert_column - Columns cannot be initialized if the coefficient field "
                                 "characteristic is not specified."));
  }
 
  static_assert(
      !isNonBasic,
      "'insert_column' not available for the chosen options. The input has to be in the form of a cell boundary.");
  matrix_.insert_column(column);
}
 
template <class PersistenceMatrixOptions>
template <class Container>
inline void Matrix<PersistenceMatrixOptions>::insert_column(const Container& column, Index columnIndex)
{
  if constexpr (!PersistenceMatrixOptions::is_z2){
    GUDHI_CHECK(colSettings_->operators.get_characteristic() != get_null_value<Characteristic>(),
                std::logic_error("Matrix::insert_column - Columns cannot be initialized if the coefficient field "
                                 "characteristic is not specified."));
  }
 
  static_assert(!isNonBasic && !PersistenceMatrixOptions::has_column_compression,
                "'insert_column' with those parameters is not available for the chosen options.");
  static_assert(!PersistenceMatrixOptions::has_row_access,
                "Columns have to be inserted at the end of the matrix when row access is enabled.");
  matrix_.insert_column(column, columnIndex);
}
 
template <class PersistenceMatrixOptions>
template <class Boundary_range>
inline typename Matrix<PersistenceMatrixOptions>::Insertion_return
Matrix<PersistenceMatrixOptions>::insert_boundary(const Boundary_range& boundary, Dimension dim)
{
  if constexpr (!PersistenceMatrixOptions::is_z2){
    GUDHI_CHECK(colSettings_->operators.get_characteristic() != get_null_value<Characteristic>(),
                std::logic_error("Matrix::insert_boundary - Columns cannot be initialized if the coefficient field "
                                 "characteristic is not specified."));
  }
 
  if constexpr (isNonBasic && !PersistenceMatrixOptions::is_of_boundary_type &&
                PersistenceMatrixOptions::column_indexation_type == Column_indexation_types::CONTAINER)
    return matrix_.insert_boundary(boundary, dim);
  else
    matrix_.insert_boundary(boundary, dim);
}
 
template <class PersistenceMatrixOptions>
template <class Boundary_range>
inline typename Matrix<PersistenceMatrixOptions>::Insertion_return
Matrix<PersistenceMatrixOptions>::insert_boundary(ID_index cellIndex,
                                                  const Boundary_range& boundary,
                                                  Dimension dim)
{
  if constexpr (!PersistenceMatrixOptions::is_z2){
    GUDHI_CHECK(colSettings_->operators.get_characteristic() != get_null_value<Characteristic>(),
                std::logic_error("Matrix::insert_boundary - Columns cannot be initialized if the coefficient field "
                                 "characteristic is not specified."));
  }
 
  static_assert(isNonBasic, "Only enabled for non-basic matrices.");
  if constexpr (!PersistenceMatrixOptions::is_of_boundary_type &&
                PersistenceMatrixOptions::column_indexation_type == Column_indexation_types::CONTAINER)
    return matrix_.insert_boundary(cellIndex, boundary, dim);
  else
    matrix_.insert_boundary(cellIndex, boundary, dim);
}
 
template <class PersistenceMatrixOptions>
inline typename Matrix<PersistenceMatrixOptions>::Returned_column& Matrix<PersistenceMatrixOptions>::get_column(
    Index columnIndex)
{
  return matrix_.get_column(columnIndex);
}
 
template <class PersistenceMatrixOptions>
inline const typename Matrix<PersistenceMatrixOptions>::Column& Matrix<PersistenceMatrixOptions>::get_column(
    Index columnIndex) const
{
  return matrix_.get_column(columnIndex);
}
 
template <class PersistenceMatrixOptions>
inline const typename Matrix<PersistenceMatrixOptions>::Column& Matrix<PersistenceMatrixOptions>::get_column(
    Index columnIndex, bool inR)
{
  // TODO: I don't think there is a particular reason why the indexation is forced, should be removed.
  static_assert(
      isNonBasic && PersistenceMatrixOptions::is_of_boundary_type &&
          (PersistenceMatrixOptions::has_vine_update || PersistenceMatrixOptions::can_retrieve_representative_cycles) &&
          PersistenceMatrixOptions::column_indexation_type != Column_indexation_types::IDENTIFIER,
      "Only enabled for position indexed RU matrices.");
 
  return matrix_.get_column(columnIndex, inR);
}
 
template <class PersistenceMatrixOptions>
inline typename Matrix<PersistenceMatrixOptions>::Returned_row& Matrix<PersistenceMatrixOptions>::get_row(
    ID_index rowIndex)
{
  static_assert(PersistenceMatrixOptions::has_row_access, "'get_row' is not available for the chosen options.");
 
  return matrix_.get_row(rowIndex);
}
 
template <class PersistenceMatrixOptions>
inline const typename Matrix<PersistenceMatrixOptions>::Row& Matrix<PersistenceMatrixOptions>::get_row(
    ID_index rowIndex) const
{
  static_assert(PersistenceMatrixOptions::has_row_access, "'get_row' is not available for the chosen options.");
 
  return matrix_.get_row(rowIndex);
}
 
template <class PersistenceMatrixOptions>
inline const typename Matrix<PersistenceMatrixOptions>::Row& Matrix<PersistenceMatrixOptions>::get_row(
    ID_index rowIndex, bool inR)
{
  static_assert(PersistenceMatrixOptions::has_row_access, "'get_row' is not available for the chosen options.");
  // TODO: I don't think there is a particular reason why the indexation is forced, should be removed.
  static_assert(
      isNonBasic && PersistenceMatrixOptions::is_of_boundary_type &&
          (PersistenceMatrixOptions::has_vine_update || PersistenceMatrixOptions::can_retrieve_representative_cycles) &&
          PersistenceMatrixOptions::column_indexation_type != Column_indexation_types::IDENTIFIER,
      "Only enabled for position indexed RU matrices.");
 
  return matrix_.get_row(rowIndex, inR);
}
 
template <class PersistenceMatrixOptions>
inline void Matrix<PersistenceMatrixOptions>::remove_column(Index columnIndex)
{
  static_assert(PersistenceMatrixOptions::has_map_column_container && !isNonBasic &&
                    !PersistenceMatrixOptions::has_column_compression,
                "'remove_column' is not available for the chosen options.");
 
  matrix_.remove_column(columnIndex);
}
 
template <class PersistenceMatrixOptions>
inline void Matrix<PersistenceMatrixOptions>::erase_empty_row(ID_index rowIndex)
{
  static_assert(
      !isNonBasic || PersistenceMatrixOptions::is_of_boundary_type || PersistenceMatrixOptions::has_removable_rows,
      "'erase_empty_row' is not available for the chosen options.");
 
  matrix_.erase_empty_row(rowIndex);
}
 
template <class PersistenceMatrixOptions>
inline void Matrix<PersistenceMatrixOptions>::remove_maximal_cell(Index columnIndex)
{
  static_assert(PersistenceMatrixOptions::has_removable_columns,
                "'remove_maximal_cell(ID_index)' is not available for the chosen options.");
  static_assert(isNonBasic && PersistenceMatrixOptions::has_vine_update,
                "'remove_maximal_cell(ID_index)' is not available for the chosen options.");
  static_assert(PersistenceMatrixOptions::is_of_boundary_type || (PersistenceMatrixOptions::has_map_column_container &&
                                                                  PersistenceMatrixOptions::has_column_pairings),
                "'remove_maximal_cell(ID_index)' is not available for the chosen options.");
 
  matrix_.remove_maximal_cell(columnIndex);
}
 
template <class PersistenceMatrixOptions>
inline void Matrix<PersistenceMatrixOptions>::remove_maximal_cell(ID_index cellIndex,
                                                                  const std::vector<ID_index>& columnsToSwap)
{
  static_assert(PersistenceMatrixOptions::has_removable_columns,
                "'remove_maximal_cell(ID_index,const std::vector<Index>&)' is not available for the chosen options.");
  static_assert(isNonBasic && !PersistenceMatrixOptions::is_of_boundary_type,
                "'remove_maximal_cell(ID_index,const std::vector<Index>&)' is not available for the chosen options.");
  static_assert(PersistenceMatrixOptions::has_map_column_container && PersistenceMatrixOptions::has_vine_update,
                "'remove_maximal_cell(ID_index,const std::vector<Index>&)' is not available for the chosen options.");
 
  matrix_.remove_maximal_cell(cellIndex, columnsToSwap);
}
 
template <class PersistenceMatrixOptions>
inline void Matrix<PersistenceMatrixOptions>::remove_last()
{
  static_assert(PersistenceMatrixOptions::has_removable_columns || !isNonBasic,
                "'remove_last' is not available for the chosen options.");
  static_assert(!PersistenceMatrixOptions::has_column_compression || isNonBasic,
                "'remove_last' is not available for the chosen options.");
  static_assert(!isNonBasic || PersistenceMatrixOptions::is_of_boundary_type ||
                    PersistenceMatrixOptions::has_map_column_container || !PersistenceMatrixOptions::has_vine_update,
                "'remove_last' is not available for the chosen options.");
 
  matrix_.remove_last();
}
 
template <class PersistenceMatrixOptions>
inline typename Matrix<PersistenceMatrixOptions>::Dimension Matrix<PersistenceMatrixOptions>::get_max_dimension()
    const
{
  static_assert(isNonBasic, "'get_max_dimension' is not available for the chosen options.");
 
  return matrix_.get_max_dimension();
}
 
template <class PersistenceMatrixOptions>
inline typename Matrix<PersistenceMatrixOptions>::Index Matrix<PersistenceMatrixOptions>::get_number_of_columns() const
{
  return matrix_.get_number_of_columns();
}
 
template <class PersistenceMatrixOptions>
inline typename Matrix<PersistenceMatrixOptions>::Dimension Matrix<PersistenceMatrixOptions>::get_column_dimension(
    Index columnIndex) const
{
  static_assert(isNonBasic, "'get_column_dimension' is not available for the chosen options.");
 
  return matrix_.get_column_dimension(columnIndex);
}
 
template <class PersistenceMatrixOptions>
template <typename Integer_index>
inline std::enable_if_t<std::is_integral_v<Integer_index> > Matrix<PersistenceMatrixOptions>::add_to(
    Integer_index sourceColumnIndex, Integer_index targetColumnIndex)
{
  matrix_.add_to(sourceColumnIndex, targetColumnIndex);
}
 
template <class PersistenceMatrixOptions>
template <class Entry_range>
inline std::enable_if_t<!std::is_integral_v<Entry_range> > Matrix<PersistenceMatrixOptions>::add_to(
    const Entry_range& sourceColumn, Index targetColumnIndex)
{
  static_assert(!isNonBasic,
                "For boundary or chain matrices, only additions with columns inside the matrix is allowed to maintain "
                "algebraic consistency.");
 
  matrix_.add_to(sourceColumn, targetColumnIndex);
}
 
template <class PersistenceMatrixOptions>
template <typename Integer_index>
inline std::enable_if_t<std::is_integral_v<Integer_index> > Matrix<PersistenceMatrixOptions>::multiply_target_and_add_to(
    Integer_index sourceColumnIndex, int coefficient, Integer_index targetColumnIndex)
{
  if constexpr (PersistenceMatrixOptions::is_z2) {
    // coef will be converted to bool, because of Element
    matrix_.multiply_target_and_add_to(sourceColumnIndex, coefficient % 2, targetColumnIndex);
  } else {
    matrix_.multiply_target_and_add_to(sourceColumnIndex, colSettings_->operators.get_value(coefficient), targetColumnIndex);
  }
}
 
template <class PersistenceMatrixOptions>
template <class Entry_range>
inline std::enable_if_t<!std::is_integral_v<Entry_range> > Matrix<PersistenceMatrixOptions>::multiply_target_and_add_to(
    const Entry_range& sourceColumn, int coefficient, Index targetColumnIndex)
{
  static_assert(!isNonBasic,
                "For boundary or chain matrices, only additions with columns inside the matrix is allowed to maintain "
                "algebraic consistency.");
 
  if constexpr (PersistenceMatrixOptions::is_z2) {
    // coef will be converted to bool, because of Element
    matrix_.multiply_target_and_add_to(sourceColumn, coefficient % 2, targetColumnIndex);
  } else {
    matrix_.multiply_target_and_add_to(sourceColumn, colSettings_->operators.get_value(coefficient), targetColumnIndex);
  }
}
 
template <class PersistenceMatrixOptions>
template <typename Integer_index>
inline std::enable_if_t<std::is_integral_v<Integer_index> > Matrix<PersistenceMatrixOptions>::multiply_source_and_add_to(
    int coefficient, Integer_index sourceColumnIndex, Integer_index targetColumnIndex)
{
  if constexpr (PersistenceMatrixOptions::is_z2) {
    // coef will be converted to bool, because of Element
    matrix_.multiply_source_and_add_to(coefficient % 2, sourceColumnIndex, targetColumnIndex);
  } else {
    matrix_.multiply_source_and_add_to(colSettings_->operators.get_value(coefficient), sourceColumnIndex, targetColumnIndex);
  }
}
 
template <class PersistenceMatrixOptions>
template <class Entry_range>
inline std::enable_if_t<!std::is_integral_v<Entry_range> > Matrix<PersistenceMatrixOptions>::multiply_source_and_add_to(
    int coefficient, const Entry_range& sourceColumn, Index targetColumnIndex)
{
  static_assert(!isNonBasic,
                "For boundary or chain matrices, only additions with columns inside the matrix is allowed to maintain "
                "algebraic consistency.");
 
  if constexpr (PersistenceMatrixOptions::is_z2) {
    // coef will be converted to bool, because of Element
    matrix_.multiply_source_and_add_to(coefficient % 2, sourceColumn, targetColumnIndex);
  } else {
    matrix_.multiply_source_and_add_to(colSettings_->operators.get_value(coefficient), sourceColumn, targetColumnIndex);
  }
}
 
template <class PersistenceMatrixOptions>
inline void Matrix<PersistenceMatrixOptions>::zero_entry(Index columnIndex, ID_index rowIndex)
{
  static_assert(PersistenceMatrixOptions::is_of_boundary_type && !PersistenceMatrixOptions::has_column_compression,
                "'zero_entry' is not available for the chosen options.");
 
  return matrix_.zero_entry(columnIndex, rowIndex);
}
 
template <class PersistenceMatrixOptions>
inline void Matrix<PersistenceMatrixOptions>::zero_entry(Index columnIndex, ID_index rowIndex, bool inR)
{
  // TODO: I don't think there is a particular reason why the indexation is forced, should be removed.
  static_assert(
      isNonBasic && PersistenceMatrixOptions::is_of_boundary_type &&
          (PersistenceMatrixOptions::has_vine_update || PersistenceMatrixOptions::can_retrieve_representative_cycles) &&
          PersistenceMatrixOptions::column_indexation_type != Column_indexation_types::IDENTIFIER,
      "Only enabled for RU matrices.");
 
  return matrix_.zero_entry(columnIndex, rowIndex, inR);
}
 
template <class PersistenceMatrixOptions>
inline void Matrix<PersistenceMatrixOptions>::zero_column(Index columnIndex)
{
  static_assert(PersistenceMatrixOptions::is_of_boundary_type && !PersistenceMatrixOptions::has_column_compression,
                "'zero_column' is not available for the chosen options.");
 
  return matrix_.zero_column(columnIndex);
}
 
template <class PersistenceMatrixOptions>
inline void Matrix<PersistenceMatrixOptions>::zero_column(Index columnIndex, bool inR)
{
  // TODO: I don't think there is a particular reason why the indexation is forced, should be removed.
  static_assert(
      isNonBasic && PersistenceMatrixOptions::is_of_boundary_type &&
          (PersistenceMatrixOptions::has_vine_update || PersistenceMatrixOptions::can_retrieve_representative_cycles) &&
          PersistenceMatrixOptions::column_indexation_type != Column_indexation_types::IDENTIFIER,
      "Only enabled for RU matrices.");
 
  return matrix_.zero_column(columnIndex, inR);
}
 
template <class PersistenceMatrixOptions>
inline bool Matrix<PersistenceMatrixOptions>::is_zero_entry(Index columnIndex, ID_index rowIndex)
{
  return matrix_.is_zero_entry(columnIndex, rowIndex);
}
 
template <class PersistenceMatrixOptions>
inline bool Matrix<PersistenceMatrixOptions>::is_zero_entry(Index columnIndex, ID_index rowIndex, bool inR) const
{
  // TODO: I don't think there is a particular reason why the indexation is forced, should be removed.
  static_assert(
      isNonBasic && PersistenceMatrixOptions::is_of_boundary_type &&
          (PersistenceMatrixOptions::has_vine_update || PersistenceMatrixOptions::can_retrieve_representative_cycles) &&
          PersistenceMatrixOptions::column_indexation_type != Column_indexation_types::IDENTIFIER,
      "Only enabled for RU matrices.");
 
  return matrix_.is_zero_entry(columnIndex, rowIndex, inR);
}
 
template <class PersistenceMatrixOptions>
inline bool Matrix<PersistenceMatrixOptions>::is_zero_column(Index columnIndex)
{
  return matrix_.is_zero_column(columnIndex);
}
 
template <class PersistenceMatrixOptions>
inline bool Matrix<PersistenceMatrixOptions>::is_zero_column(Index columnIndex, bool inR)
{
  // TODO: I don't think there is a particular reason why the indexation is forced, should be removed.
  static_assert(
      isNonBasic && PersistenceMatrixOptions::is_of_boundary_type &&
          (PersistenceMatrixOptions::has_vine_update || PersistenceMatrixOptions::can_retrieve_representative_cycles) &&
          PersistenceMatrixOptions::column_indexation_type != Column_indexation_types::IDENTIFIER,
      "Only enabled for RU matrices.");
 
  return matrix_.is_zero_column(columnIndex, inR);
}
 
template <class PersistenceMatrixOptions>
inline typename Matrix<PersistenceMatrixOptions>::Index Matrix<PersistenceMatrixOptions>::get_column_with_pivot(
    ID_index cellIndex) const
{
  static_assert(isNonBasic && (!PersistenceMatrixOptions::is_of_boundary_type ||
                               (PersistenceMatrixOptions::has_vine_update ||
                                PersistenceMatrixOptions::can_retrieve_representative_cycles)),
                "'get_column_with_pivot' is not available for the chosen options.");
 
  return matrix_.get_column_with_pivot(cellIndex);
}
 
template <class PersistenceMatrixOptions>
inline typename Matrix<PersistenceMatrixOptions>::ID_index Matrix<PersistenceMatrixOptions>::get_pivot(
    Index columnIndex)
{
  static_assert(isNonBasic, "'get_pivot' is not available for the chosen options.");
 
  return matrix_.get_pivot(columnIndex);
}
 
template <class PersistenceMatrixOptions>
inline Matrix<PersistenceMatrixOptions>& Matrix<PersistenceMatrixOptions>::operator=(Matrix other)
{
  swap(matrix_, other.matrix_);
  std::swap(colSettings_, other.colSettings_);
 
  return *this;
}
 
template <class PersistenceMatrixOptions>
inline void Matrix<PersistenceMatrixOptions>::print()
{
  return matrix_.print();
}
 
template <class PersistenceMatrixOptions>
inline const typename Matrix<PersistenceMatrixOptions>::Barcode&
Matrix<PersistenceMatrixOptions>::get_current_barcode()
{
  static_assert(PersistenceMatrixOptions::has_column_pairings, "This method was not enabled.");
 
  return matrix_.get_current_barcode();
}
 
template <class PersistenceMatrixOptions>
inline const typename Matrix<PersistenceMatrixOptions>::Barcode&
Matrix<PersistenceMatrixOptions>::get_current_barcode() const
{
  static_assert(PersistenceMatrixOptions::has_column_pairings, "This method was not enabled.");
  static_assert(
      !PersistenceMatrixOptions::is_of_boundary_type || PersistenceMatrixOptions::has_vine_update ||
          PersistenceMatrixOptions::can_retrieve_representative_cycles,
      "'get_current_barcode' is not const for boundary matrices as the barcode is only computed when explicitly "
      "asked.");
 
  return matrix_.get_current_barcode();
}
 
template <class PersistenceMatrixOptions>
inline void Matrix<PersistenceMatrixOptions>::swap_columns(Index columnIndex1, Index columnIndex2)
{
  static_assert(
      (!isNonBasic && !PersistenceMatrixOptions::has_column_compression) ||
          (isNonBasic && PersistenceMatrixOptions::is_of_boundary_type && !PersistenceMatrixOptions::has_vine_update &&
           !PersistenceMatrixOptions::can_retrieve_representative_cycles),
      "This method was not enabled.");
  return matrix_.swap_columns(columnIndex1, columnIndex2);
}
 
template <class PersistenceMatrixOptions>
inline void Matrix<PersistenceMatrixOptions>::swap_rows(Index rowIndex1, Index rowIndex2)
{
  static_assert(
      (!isNonBasic && !PersistenceMatrixOptions::has_column_compression) ||
          (isNonBasic && PersistenceMatrixOptions::is_of_boundary_type && !PersistenceMatrixOptions::has_vine_update &&
           !PersistenceMatrixOptions::can_retrieve_representative_cycles),
      "This method was not enabled.");
  return matrix_.swap_rows(rowIndex1, rowIndex2);
}
 
template <class PersistenceMatrixOptions>
inline bool Matrix<PersistenceMatrixOptions>::vine_swap_with_z_eq_1_case(Pos_index index)
{
  static_assert(PersistenceMatrixOptions::has_vine_update, "This method was not enabled.");
  static_assert(PersistenceMatrixOptions::column_indexation_type == Column_indexation_types::POSITION ||
                    (PersistenceMatrixOptions::is_of_boundary_type &&
                     PersistenceMatrixOptions::column_indexation_type == Column_indexation_types::CONTAINER),
                "This method was not enabled.");
  return matrix_.vine_swap_with_z_eq_1_case(index);
}
 
template <class PersistenceMatrixOptions>
inline typename Matrix<PersistenceMatrixOptions>::Index Matrix<PersistenceMatrixOptions>::vine_swap_with_z_eq_1_case(
    Index columnIndex1, Index columnIndex2)
{
  static_assert(PersistenceMatrixOptions::has_vine_update, "This method was not enabled.");
  static_assert(PersistenceMatrixOptions::column_indexation_type == Column_indexation_types::IDENTIFIER ||
                    (!PersistenceMatrixOptions::is_of_boundary_type &&
                     PersistenceMatrixOptions::column_indexation_type == Column_indexation_types::CONTAINER),
                "This method was not enabled.");
 
  return matrix_.vine_swap_with_z_eq_1_case(columnIndex1, columnIndex2);
}
 
template <class PersistenceMatrixOptions>
inline bool Matrix<PersistenceMatrixOptions>::vine_swap(Pos_index index)
{
  static_assert(PersistenceMatrixOptions::has_vine_update, "This method was not enabled.");
  static_assert(PersistenceMatrixOptions::column_indexation_type == Column_indexation_types::POSITION ||
                    (PersistenceMatrixOptions::is_of_boundary_type &&
                     PersistenceMatrixOptions::column_indexation_type == Column_indexation_types::CONTAINER),
                "This method was not enabled.");
  return matrix_.vine_swap(index);
}
 
template <class PersistenceMatrixOptions>
inline typename Matrix<PersistenceMatrixOptions>::Index Matrix<PersistenceMatrixOptions>::vine_swap(
    Index columnIndex1, Index columnIndex2)
{
  static_assert(PersistenceMatrixOptions::has_vine_update, "This method was not enabled.");
  static_assert(PersistenceMatrixOptions::column_indexation_type == Column_indexation_types::IDENTIFIER ||
                    (!PersistenceMatrixOptions::is_of_boundary_type &&
                     PersistenceMatrixOptions::column_indexation_type == Column_indexation_types::CONTAINER),
                "This method was not enabled.");
  return matrix_.vine_swap(columnIndex1, columnIndex2);
}
 
template <class PersistenceMatrixOptions>
inline void Matrix<PersistenceMatrixOptions>::update_representative_cycles()
{
  static_assert(PersistenceMatrixOptions::can_retrieve_representative_cycles, "This method was not enabled.");
  matrix_.update_representative_cycles();
}
 
template <class PersistenceMatrixOptions>
inline const std::vector<typename Matrix<PersistenceMatrixOptions>::Cycle>&
Matrix<PersistenceMatrixOptions>::get_representative_cycles()
{
  static_assert(PersistenceMatrixOptions::can_retrieve_representative_cycles, "This method was not enabled.");
  return matrix_.get_representative_cycles();
}
 
template <class PersistenceMatrixOptions>
inline const typename Matrix<PersistenceMatrixOptions>::Cycle&
Matrix<PersistenceMatrixOptions>::get_representative_cycle(const Bar& bar)
{
  static_assert(PersistenceMatrixOptions::can_retrieve_representative_cycles, "This method was not enabled.");
  return matrix_.get_representative_cycle(bar);
}
 
template <class PersistenceMatrixOptions>
inline constexpr void Matrix<PersistenceMatrixOptions>::_assert_options()
{
  static_assert(
      PersistenceMatrixOptions::column_type != Column_types::HEAP || !PersistenceMatrixOptions::has_row_access,
      "Row access is not possible for heap columns.");
  static_assert(!PersistenceMatrixOptions::has_vine_update || PersistenceMatrixOptions::is_z2,
                "Vine update currently works only for Z_2 coefficients.");
  // static_assert(!PersistenceMatrixOptions::can_retrieve_representative_cycles || PersistenceMatrixOptions::is_z2,
  //               "Representative cycles can currently only be computed with Z_2 coefficients.");
  static_assert(
      PersistenceMatrixOptions::column_type != Column_types::HEAP || !PersistenceMatrixOptions::has_column_compression,
      "Column compression not compatible with heap columns.");
 
  // // This should be warnings instead, as PersistenceMatrixOptions::has_column_compression is just ignored in those
  // // cases and don't produces errors as long as the corresponding methods are not called.
  // static_assert(!PersistenceMatrixOptions::has_column_compression || !PersistenceMatrixOptions::has_column_pairings,
  //               "Column compression not available to compute persistence homology (it would bring no advantages; "
  //               "use it for co-homology instead).");
  // static_assert(!PersistenceMatrixOptions::has_column_compression || !PersistenceMatrixOptions::has_vine_update,
  //               "Column compression not available for vineyards.");
  // static_assert(!PersistenceMatrixOptions::has_column_compression ||
  //                   !PersistenceMatrixOptions::can_retrieve_representative_cycles,
  //               "Column compression not available to retrieve representative cycles.");
  // // Would column removal while column compression be useful? If yes, should erase() remove a single column or the
  // // class of columns identical to the input?
  // // For a single column, I have an implementation for union-find (not the current one) which allows deleting a
  // // single element in constant time, but find becomes log n in worst case.
  // // For a column class, we either just ignore the removed class (constant time), but this results in memory
  // // residues, or, we have an erase method which is at least linear in the size of the class.
  // static_assert(
  //     !PersistenceMatrixOptions::has_column_compression || !PersistenceMatrixOptions::has_map_column_container,
  //     "When column compression is used, the removal of columns is not implemented yet.");
}
 
}  // namespace persistence_matrix
}  // namespace Gudhi
 
#endif  // MASTER_MATRIX_H