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/Fields/Z2_field_operators.h>
 
#include <gudhi/Persistence_matrix/overlay_ididx_to_matidx.h>
#include <gudhi/Persistence_matrix/overlay_posidx_to_matidx.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/cell_constructors.h>
#include <gudhi/Persistence_matrix/columns/cell_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_type;                 
  using id_index = typename PersistenceMatrixOptions::index_type;              
  using pos_index = typename PersistenceMatrixOptions::index_type;             
  using dimension_type = typename PersistenceMatrixOptions::dimension_type;    
  using Field_operators =
      typename std::conditional<PersistenceMatrixOptions::is_z2, 
                                Gudhi::persistence_fields::Z2_field_operators, 
                                typename PersistenceMatrixOptions::Field_coeff_operators
                               >::type;
  using element_type = typename Field_operators::element_type;
  using characteristic_type = typename Field_operators::characteristic_type;
 
  // TODO: move outside? unify with other bar types in Gudhi?
  struct Bar {
    Bar() : dim(-1), birth(-1), death(-1) {}
 
    Bar(dimension_type dim, pos_index birth, pos_index death) : dim(dim), birth(birth), death(death) {}
 
    dimension_type dim; 
    pos_index birth;    
    pos_index death;    
    inline friend std::ostream &operator<<(std::ostream &stream, const Bar& bar) {
      stream << "[" << bar.dim << "] ";
      stream << bar.birth << ", " << bar.death;
      stream << "\n";
      return stream;
    }
  };
 
  //tags for boost to associate a row and a column to a same cell
  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 a cell can inherit both a row and a column hook.
  struct Dummy_row_hook {};
  struct Dummy_column_hook {};
 
  using row_hook_type = typename std::conditional<PersistenceMatrixOptions::has_row_access &&
                                                      PersistenceMatrixOptions::has_intrusive_rows,
                                                  base_hook_matrix_row, 
                                                  Dummy_row_hook
                                                 >::type;
  using column_hook_type = 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 cell (additionnaly to the row index). Necessary only with row access.
  using Cell_column_index_option =
      typename std::conditional<PersistenceMatrixOptions::has_row_access,
                                Cell_column_index<index>,
                                Dummy_cell_column_index_mixin
                               >::type;
  //Option to store the value of the cell. 
  //Unnecessary for values in Z_2 as there are always 1 (0-valued cells are never stored).
  using Cell_field_element_option =
      typename std::conditional<PersistenceMatrixOptions::is_z2,
                                Dummy_cell_field_element_mixin,
                                Cell_field_element<element_type>
                               >::type;
  using Cell_type = Cell<Matrix<PersistenceMatrixOptions> >;
 
  inline static New_cell_constructor<Cell_type> defaultCellConstructor;
  using Cell_constructor = Pool_cell_constructor<Cell_type>;
 
  using cell_rep_type = typename std::conditional<PersistenceMatrixOptions::is_z2,
                                                  id_index,
                                                  std::pair<id_index, element_type>
                                                 >::type;
 
  struct RowCellComp {
    bool operator()(const Cell_type& c1, const Cell_type& c2) const {
      return c1.get_column_index() < c2.get_column_index();
    }
  };
 
  using Row_type =
      typename std::conditional<PersistenceMatrixOptions::has_intrusive_rows,
                                boost::intrusive::list<Cell_type, 
                                                       boost::intrusive::constant_time_size<false>,
                                                       boost::intrusive::base_hook<base_hook_matrix_row>
                                                      >,
                                std::set<Cell_type, RowCellComp>
                               >::type;
 
  using row_container_type =
      typename std::conditional<PersistenceMatrixOptions::has_removable_rows,
                                std::map<id_index, Row_type>,
                                std::vector<Row_type>
                               >::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_type, 
                                                  row_container_type, 
                                                  PersistenceMatrixOptions::has_removable_rows, 
                                                  id_index>,
                                Dummy_matrix_row_access
                               >::type;
 
  template <typename value_type>
  using dictionnary_type =
      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 Heap_column_type = Heap_column<Matrix<PersistenceMatrixOptions> >;
  using List_column_type = List_column<Matrix<PersistenceMatrixOptions> >;
  using Vector_column_type = Vector_column<Matrix<PersistenceMatrixOptions> >;
  using Naive_vector_column_type = Naive_vector_column<Matrix<PersistenceMatrixOptions> >;
  using Set_column_type = Set_column<Matrix<PersistenceMatrixOptions> >;
  using Unordered_set_column_type = Unordered_set_column<Matrix<PersistenceMatrixOptions> >;
  using Intrusive_list_column_type = Intrusive_list_column<Matrix<PersistenceMatrixOptions> >;
  using Intrusive_set_column_type = Intrusive_set_column<Matrix<PersistenceMatrixOptions> >;
 
  using Column_type = typename std::conditional<
        PersistenceMatrixOptions::column_type == Column_types::HEAP, 
        Heap_column_type,
        typename std::conditional<
            PersistenceMatrixOptions::column_type == Column_types::LIST, 
            List_column_type,
            typename std::conditional<
                PersistenceMatrixOptions::column_type == Column_types::SET, 
                Set_column_type,
                typename std::conditional<
                    PersistenceMatrixOptions::column_type == Column_types::UNORDERED_SET, 
                    Unordered_set_column_type,
                    typename std::conditional<
                        PersistenceMatrixOptions::column_type == Column_types::VECTOR, 
                        Vector_column_type,
                        typename std::conditional<
                            PersistenceMatrixOptions::column_type == Column_types::INTRUSIVE_LIST, 
                            Intrusive_list_column_type,
                            typename std::conditional<
                                PersistenceMatrixOptions::column_type == Column_types::NAIVE_VECTOR,
                                Naive_vector_column_type, 
                                Intrusive_set_column_type
                                >::type
                            >::type
                        >::type
                    >::type
                >::type
            >::type
        >::type;
 
  struct Column_z2_settings{
    Column_z2_settings() : cellConstructor() {}
    Column_z2_settings([[maybe_unused]] characteristic_type characteristic) : cellConstructor() {}
    Column_z2_settings(const Column_z2_settings& toCopy) : cellConstructor() {}
 
    Cell_constructor cellConstructor;   //will be replaced by more specific allocators depending on the column type.
  };
 
  struct Column_zp_settings {
    Column_zp_settings() : operators(), cellConstructor() {}
    //purposely triggers operators() instead of operators(characteristic) as the "dummy" values for the different
    //operators can be different from -1.
    Column_zp_settings(characteristic_type characteristic) : operators(), cellConstructor() {
      if (characteristic != static_cast<characteristic_type>(-1)) operators.set_characteristic(characteristic);
    }
    Column_zp_settings(const Column_zp_settings& toCopy)
        : operators(toCopy.operators.get_characteristic()), cellConstructor() {}
 
    Field_operators operators;
    Cell_constructor cellConstructor;   //will be replaced by more specific allocators depending on the column type.
  };
 
  // struct Column_z2_with_rows_settings {
  //   Column_z2_with_rows_settings() : cellConstructor(), rows(nullptr) {}
  //   Column_z2_with_rows_settings([[maybe_unused]] characteristic_type characteristic)
  //       : cellConstructor(), rows(nullptr) {}
 
  //   Cell_constructor cellConstructor;
  //   row_container_type* rows;
  // };
 
  // struct Column_zp_with_rows_settings {
  //   Column_zp_with_rows_settings() : operators(), cellConstructor(), rows(nullptr) {}
  //   Column_zp_with_rows_settings(characteristic_type characteristic)
  //       : operators(characteristic), cellConstructor(), rows(nullptr) {}
 
  //   Field_operators operators;
  //   Cell_constructor cellConstructor;
  //   row_container_type* 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_type =
      typename std::conditional<PersistenceMatrixOptions::has_map_column_container, 
                                std::unordered_map<index, Column_type>,
                                std::vector<Column_type>
                               >::type;
 
  static const bool hasFixedBarcode = Option_list::is_of_boundary_type && !PersistenceMatrixOptions::has_vine_update;
  using barcode_type =
      typename std::conditional<hasFixedBarcode, 
                                std::vector<Bar>, 
                                typename std::conditional<PersistenceMatrixOptions::has_removable_columns, 
                                  std::list<Bar>, 
                                  std::vector<Bar>
                                >::type
                               >::type;
  using bar_dictionnary_type = 
      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_type::iterator>,
                                  std::vector<index>
                                >::type
                               >::type;
 
  //default type for boundaries to permit list initialization directly in function parameters
  using boundary_type = typename std::conditional<PersistenceMatrixOptions::is_z2, 
                                                  std::initializer_list<id_index>,
                                                  std::initializer_list<std::pair<id_index, element_type> >
                                                 >::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_type>,
                                Matrix_max_dimension_holder<dimension_type>
                               >::type,
      Dummy_matrix_dimension_holder
    >::type;
 
  using Base_matrix_type =
      typename std::conditional<PersistenceMatrixOptions::has_column_compression, 
                                Base_matrix_with_column_compression<Matrix<PersistenceMatrixOptions> >,
                                Base_matrix<Matrix<PersistenceMatrixOptions> >
                               >::type;
  using Boundary_matrix_type = Boundary_matrix<Matrix<PersistenceMatrixOptions> >;
  using RU_matrix_type = RU_matrix<Matrix<PersistenceMatrixOptions> >;
  using Chain_matrix_type = 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_type = 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_type =
      typename std::conditional<!isNonBasic && PersistenceMatrixOptions::has_column_compression,
                                const Column_type,
                                Column_type
                               >::type;
  //The returned row is `const` if the matrix uses column compression
  using returned_row_type =
      typename std::conditional<!isNonBasic && PersistenceMatrixOptions::has_column_compression,
                                const Row_type,
                                Row_type
                               >::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_type =
      typename std::conditional<PersistenceMatrixOptions::is_of_boundary_type || !isNonBasic ||
                                    PersistenceMatrixOptions::column_indexation_type ==
                                        Column_indexation_types::POSITION,
                                void, 
                                std::vector<cell_rep_type>
                               >::type;
 
  Matrix();
  template <class Container_type = boundary_type>
  Matrix(const std::vector<Container_type>& columns, characteristic_type characteristic = 11);
  Matrix(int numberOfColumns, characteristic_type characteristic = static_cast<characteristic_type>(-1));
  Matrix(const std::function<bool(pos_index,pos_index)>& birthComparator, 
         const std::function<bool(pos_index,pos_index)>& deathComparator);
  template <class Boundary_type = boundary_type>
  Matrix(const std::vector<Boundary_type>& orderedBoundaries, 
         const std::function<bool(pos_index,pos_index)>& birthComparator,
         const std::function<bool(pos_index,pos_index)>& deathComparator, 
         characteristic_type 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_type characteristic = static_cast<characteristic_type>(-1));
  Matrix(const Matrix& matrixToCopy);
  Matrix(Matrix&& other) noexcept;
 
  ~Matrix();
 
  //TODO: compatibily with multi fields:
  //  - set_characteristic(characteristic_type min, characteristic_type max)
  //  - readapt reduction?
  void set_characteristic(characteristic_type 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_type>
  void insert_column(const Container_type& column);
  template <class Container_type>
  void insert_column(const Container_type& 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_type = boundary_type>
  insertion_return_type insert_boundary(const Boundary_type& boundary, dimension_type dim = -1);
  template <class Boundary_type>
  insertion_return_type insert_boundary(id_index faceIndex, const Boundary_type& boundary, dimension_type dim = -1);
 
  returned_column_type& get_column(index columnIndex);
  const Column_type& 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_type& get_column(index columnIndex, bool inR);
 
  //TODO: update column indices when reordering rows (after lazy swap) such that always MatIdx are returned.
  returned_row_type& get_row(id_index rowIndex);
  const Row_type& 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_type& 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 homogenise.
  void remove_maximal_face(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 (face, {})...
  void remove_maximal_face(id_index faceIndex, const std::vector<id_index>& columnsToSwap);
  void remove_last();
 
  dimension_type get_max_dimension() const;
  index get_number_of_columns() const;
  dimension_type get_column_dimension(index columnIndex) const;
 
  template <typename Index_type>
  std::enable_if_t<std::is_integral_v<Index_type> > add_to(Index_type sourceColumnIndex, Index_type targetColumnIndex);
  template <class Cell_range>
  std::enable_if_t<!std::is_integral_v<Cell_range> > add_to(const Cell_range& sourceColumn, index targetColumnIndex);
 
  template <typename Index_type>
  std::enable_if_t<std::is_integral_v<Index_type> > multiply_target_and_add_to(Index_type sourceColumnIndex,
                                                                               int coefficient,
                                                                               Index_type targetColumnIndex);
  template <class Cell_range>
  std::enable_if_t<!std::is_integral_v<Cell_range> > multiply_target_and_add_to(const Cell_range& sourceColumn,
                                                                                int coefficient,
                                                                                index targetColumnIndex);
 
  template <typename Index_type>
  std::enable_if_t<std::is_integral_v<Index_type> > multiply_source_and_add_to(int coefficient,
                                                                               Index_type sourceColumnIndex,
                                                                               Index_type targetColumnIndex);
  template <class Cell_range>
  std::enable_if_t<!std::is_integral_v<Cell_range> > multiply_source_and_add_to(int coefficient,
                                                                                const Cell_range& sourceColumn,
                                                                                index targetColumnIndex);
 
  void zero_cell(index columnIndex, id_index rowIndex);
  void zero_cell(index columnIndex, id_index rowIndex, bool inR);
  void zero_column(index columnIndex);
  void zero_column(index columnIndex, bool inR);
  bool is_zero_cell(index columnIndex, id_index rowIndex);
  bool is_zero_cell(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 faceIndex) 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_type& get_current_barcode();
  const barcode_type& 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_type>& get_representative_cycles();
  const cycle_type& get_representative_cycle(const Bar& bar);
 
 private:
  using Matrix_type = 
    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,
                    RU_matrix_type, Id_to_index_overlay<RU_matrix_type, Matrix<PersistenceMatrixOptions> >
                >::type,
                typename std::conditional<
                    PersistenceMatrixOptions::column_indexation_type == Column_indexation_types::CONTAINER || 
                      PersistenceMatrixOptions::column_indexation_type == Column_indexation_types::POSITION,
                    Boundary_matrix_type,
                    Id_to_index_overlay<Boundary_matrix_type, Matrix<PersistenceMatrixOptions> >
                >::type
            >::type,
            typename std::conditional<
                PersistenceMatrixOptions::column_indexation_type == Column_indexation_types::CONTAINER, 
                Chain_matrix_type,
                typename std::conditional<
                    PersistenceMatrixOptions::column_indexation_type == Column_indexation_types::POSITION,
                    Position_to_index_overlay<Chain_matrix_type, Matrix<PersistenceMatrixOptions> >,
                    Id_to_index_overlay<Chain_matrix_type, Matrix<PersistenceMatrixOptions> >
                >::type
            >::type
        >::type,
        Base_matrix_type
    >::type;
 
  // Field_operators* operators_;
  // Cell_constructor* cellPool_;
  Column_settings* colSettings_;  //pointer because the of swap operator on matrix_ which also stores the pointer
  Matrix_type 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, comparaison functions for the columns have to be provided.");
  _assert_options();
}
 
template <class PersistenceMatrixOptions>
template <class Container_type>
inline Matrix<PersistenceMatrixOptions>::Matrix(const std::vector<Container_type>& columns,
                                                characteristic_type 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, comparaison functions for the columns "
                "have to be provided.");
  _assert_options();
}
 
template <class PersistenceMatrixOptions>
inline Matrix<PersistenceMatrixOptions>::Matrix(int numberOfColumns, characteristic_type 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, comparaison 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 barcodes are not recorded.");
  _assert_options();
}
 
template <class PersistenceMatrixOptions>
template <class Boundary_type>
inline Matrix<PersistenceMatrixOptions>::Matrix(const std::vector<Boundary_type>& orderedBoundaries,
                               const std::function<bool(pos_index,pos_index)>& birthComparator, 
                               const std::function<bool(pos_index,pos_index)>& deathComparator,
                               characteristic_type 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 barcodes are 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_type 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 barcodes are 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_type characteristic) 
{
  if constexpr (!PersistenceMatrixOptions::is_z2) {
    if (colSettings_->operators.get_characteristic() != static_cast<characteristic_type>(-1)) {
      std::cerr << "Warning: Characteristic already initialised. Changing it could lead to incoherences in the matrice "
                   "as the modulo was already applied to values in existing columns.";
    }
 
    colSettings_->operators.set_characteristic(characteristic);
  }
}
 
template <class PersistenceMatrixOptions>
template <class Container_type>
inline void Matrix<PersistenceMatrixOptions>::insert_column(const Container_type& column) 
{
  if constexpr (!PersistenceMatrixOptions::is_z2){
    GUDHI_CHECK(colSettings_->operators.get_characteristic() != static_cast<characteristic_type>(-1),
                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 face boundary.");
  matrix_.insert_column(column);
}
 
template <class PersistenceMatrixOptions>
template <class Container_type>
inline void Matrix<PersistenceMatrixOptions>::insert_column(const Container_type& column, index columnIndex) 
{
  if constexpr (!PersistenceMatrixOptions::is_z2){
    GUDHI_CHECK(colSettings_->operators.get_characteristic() != static_cast<characteristic_type>(-1),
                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_type>
inline typename Matrix<PersistenceMatrixOptions>::insertion_return_type
Matrix<PersistenceMatrixOptions>::insert_boundary(const Boundary_type& boundary, dimension_type dim)
{
  if constexpr (!PersistenceMatrixOptions::is_z2){
    GUDHI_CHECK(colSettings_->operators.get_characteristic() != static_cast<characteristic_type>(-1),
                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_type>
inline typename Matrix<PersistenceMatrixOptions>::insertion_return_type
Matrix<PersistenceMatrixOptions>::insert_boundary(id_index faceIndex, 
                                                  const Boundary_type& boundary,
                                                  dimension_type dim)
{
  if constexpr (!PersistenceMatrixOptions::is_z2){
    GUDHI_CHECK(colSettings_->operators.get_characteristic() != static_cast<characteristic_type>(-1),
                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(faceIndex, boundary, dim);
  else
    matrix_.insert_boundary(faceIndex, boundary, dim);
}
 
template <class PersistenceMatrixOptions>
inline typename Matrix<PersistenceMatrixOptions>::returned_column_type& Matrix<PersistenceMatrixOptions>::get_column(
    index columnIndex)
{
  return matrix_.get_column(columnIndex);
}
 
template <class PersistenceMatrixOptions>
inline const typename Matrix<PersistenceMatrixOptions>::Column_type& Matrix<PersistenceMatrixOptions>::get_column(
    index columnIndex) const 
{
  return matrix_.get_column(columnIndex);
}
 
template <class PersistenceMatrixOptions>
inline const typename Matrix<PersistenceMatrixOptions>::Column_type& 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_type& 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_type& 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_type& 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_face(index columnIndex) 
{
  static_assert(PersistenceMatrixOptions::has_removable_columns,
                "'remove_maximal_face(id_index)' is not available for the chosen options.");
  static_assert(isNonBasic && PersistenceMatrixOptions::has_vine_update,
                "'remove_maximal_face(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_face(id_index)' is not available for the chosen options.");
 
  matrix_.remove_maximal_face(columnIndex);
}
 
template <class PersistenceMatrixOptions>
inline void Matrix<PersistenceMatrixOptions>::remove_maximal_face(id_index faceIndex,
                                                                  const std::vector<id_index>& columnsToSwap)
{
  static_assert(PersistenceMatrixOptions::has_removable_columns,
                "'remove_maximal_face(id_index,const std::vector<index>&)' is not available for the chosen options.");
  static_assert(isNonBasic && !PersistenceMatrixOptions::is_of_boundary_type,
                "'remove_maximal_face(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_face(id_index,const std::vector<index>&)' is not available for the chosen options.");
 
  matrix_.remove_maximal_face(faceIndex, 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_type 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_type 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 Index_type>
inline std::enable_if_t<std::is_integral_v<Index_type> > Matrix<PersistenceMatrixOptions>::add_to(
    Index_type sourceColumnIndex, Index_type targetColumnIndex)
{
  matrix_.add_to(sourceColumnIndex, targetColumnIndex);
}
 
template <class PersistenceMatrixOptions>
template <class Cell_range>
inline std::enable_if_t<!std::is_integral_v<Cell_range> > Matrix<PersistenceMatrixOptions>::add_to(
    const Cell_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 Index_type>
inline std::enable_if_t<std::is_integral_v<Index_type> > Matrix<PersistenceMatrixOptions>::multiply_target_and_add_to(
    Index_type sourceColumnIndex, int coefficient, Index_type targetColumnIndex) 
{
  if constexpr (PersistenceMatrixOptions::is_z2) {
    // coef will be converted to bool, because of element_type
    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 Cell_range>
inline std::enable_if_t<!std::is_integral_v<Cell_range> > Matrix<PersistenceMatrixOptions>::multiply_target_and_add_to(
    const Cell_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_type
    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 Index_type>
inline std::enable_if_t<std::is_integral_v<Index_type> > Matrix<PersistenceMatrixOptions>::multiply_source_and_add_to(
    int coefficient, Index_type sourceColumnIndex, Index_type targetColumnIndex) 
{
  if constexpr (PersistenceMatrixOptions::is_z2) {
    // coef will be converted to bool, because of element_type
    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 Cell_range>
inline std::enable_if_t<!std::is_integral_v<Cell_range> > Matrix<PersistenceMatrixOptions>::multiply_source_and_add_to(
    int coefficient, const Cell_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_type
    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_cell(index columnIndex, id_index rowIndex) 
{
  static_assert(PersistenceMatrixOptions::is_of_boundary_type && !PersistenceMatrixOptions::has_column_compression,
                "'zero_cell' is not available for the chosen options.");
 
  return matrix_.zero_cell(columnIndex, rowIndex);
}
 
template <class PersistenceMatrixOptions>
inline void Matrix<PersistenceMatrixOptions>::zero_cell(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_cell(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_cell(index columnIndex, id_index rowIndex) 
{
  return matrix_.is_zero_cell(columnIndex, rowIndex);
}
 
template <class PersistenceMatrixOptions>
inline bool Matrix<PersistenceMatrixOptions>::is_zero_cell(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_cell(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 faceIndex) 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(faceIndex);
}
 
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_type&
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_type&
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 explicitely "
      "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_type>&
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_type&
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,
  //               "Representaive 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