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 Inria
 *
 *    Modification(s):
 *      - 2025/11 Jānis Lazovskis: Added insert_maximal_cell method
 *      - 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 = std::conditional_t<PersistenceMatrixOptions::is_z2,
                                             Gudhi::persistence_fields::Z2_field_operators,
                                             typename PersistenceMatrixOptions::Field_coeff_operators>;
  using Element = typename Field_operators::Element;
  using Characteristic = typename Field_operators::Characteristic;

  template <typename T>
  static constexpr 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 =
      std::conditional_t<PersistenceMatrixOptions::has_row_access && PersistenceMatrixOptions::has_intrusive_rows,
                         Base_hook_matrix_row,
                         Dummy_row_hook>;
  using Column_hook =
      std::conditional_t<PersistenceMatrixOptions::column_type == Column_types::INTRUSIVE_LIST,
                         Base_hook_matrix_list_column,
                         std::conditional_t<PersistenceMatrixOptions::column_type == Column_types::INTRUSIVE_SET,
                                            Base_hook_matrix_set_column,
                                            Dummy_column_hook> >;
 
  // Option to store the column index within the entry (additionally to the row index). Necessary only with row access.
  using Entry_column_index_option = std::conditional_t<PersistenceMatrixOptions::has_row_access,
                                                       Entry_column_index<Index>,
                                                       Dummy_entry_column_index_mixin>;
  // 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 = std::
      conditional_t<PersistenceMatrixOptions::is_z2, Dummy_entry_field_element_mixin, Entry_field_element<Element> >;
  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 =
      std::conditional_t<PersistenceMatrixOptions::is_z2, ID_index, std::pair<ID_index, Element> >;

  static ID_index get_row_index(const Matrix_entry& e)
  {
    return e.get_row_index();
  }
  static ID_index get_row_index(const Entry_representative& e)
  {
    if constexpr (PersistenceMatrixOptions::is_z2) {
      return e;
    } else {
      return e.first;
    }
  }
  static ID_index& get_row_index(Entry_representative& e)
  {
    if constexpr (PersistenceMatrixOptions::is_z2) {
      return e;
    } else {
      return e.first;
    }
  }
  static Element get_element(const Matrix_entry& e)
  {
    return e.get_element();
  }
  static Element get_element(const Entry_representative& e)
  {
    if constexpr (PersistenceMatrixOptions::is_z2) {
      return Field_operators::get_multiplicative_identity();
    } else {
      return e.second;
    }
  }

  struct RowEntryComp {
    bool operator()(const Matrix_entry& c1, const Matrix_entry& c2) const
    {
      return c1.get_column_index() < c2.get_column_index();
    }
  };

  using Row = std::conditional_t<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> >;
 
  using Row_container =
      std::conditional_t<PersistenceMatrixOptions::has_removable_rows, std::map<ID_index, Row>, std::vector<Row> >;
 
  // Row access at column level
  using Row_access_option = std::conditional_t<PersistenceMatrixOptions::has_row_access,
                                               Row_access<Matrix<PersistenceMatrixOptions> >,
                                               Dummy_row_access>;
  // Row access at matrix level
  using Matrix_row_access_option =
      std::conditional_t<PersistenceMatrixOptions::has_row_access,
                         Matrix_row_access<Row, Row_container, PersistenceMatrixOptions::has_removable_rows, ID_index>,
                         Dummy_matrix_row_access>;
 
  template <typename value_type>
  using Dictionary = std::conditional_t<PersistenceMatrixOptions::has_map_column_container,
                                        std::unordered_map<unsigned int, value_type>,
                                        std::vector<value_type> >;
 
  static const bool isNonBasic = PersistenceMatrixOptions::has_column_pairings ||
                                 PersistenceMatrixOptions::has_vine_update ||
                                 PersistenceMatrixOptions::can_retrieve_representative_cycles;
 
  using Column_dimension_option = std::
      conditional_t<isNonBasic, Column_dimension_holder<Matrix<PersistenceMatrixOptions> >, Dummy_dimension_holder>;
  // Extra information needed for a column when the matrix is a @ref chainmatrix "chain matrix".
  using Chain_column_option = std::conditional_t<isNonBasic && !PersistenceMatrixOptions::is_of_boundary_type,
                                                 Chain_column_extra_properties<Matrix<PersistenceMatrixOptions> >,
                                                 Dummy_chain_properties>;
 
  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 = std::conditional_t<
      PersistenceMatrixOptions::column_type == Column_types::HEAP,
      Matrix_heap_column,
      std::conditional_t<
          PersistenceMatrixOptions::column_type == Column_types::LIST,
          Matrix_list_column,
          std::conditional_t<
              PersistenceMatrixOptions::column_type == Column_types::SET,
              Matrix_set_column,
              std::conditional_t<
                  PersistenceMatrixOptions::column_type == Column_types::UNORDERED_SET,
                  Matrix_unordered_set_column,
                  std::conditional_t<
                      PersistenceMatrixOptions::column_type == Column_types::VECTOR,
                      Matrix_vector_column,
                      std::conditional_t<
                          PersistenceMatrixOptions::column_type == Column_types::INTRUSIVE_LIST,
                          Matrix_intrusive_list_column,
                          std::conditional_t<
                              PersistenceMatrixOptions::column_type == Column_types::NAIVE_VECTOR,
                              Matrix_naive_vector_column,
                              std::conditional_t<PersistenceMatrixOptions::column_type == Column_types::SMALL_VECTOR,
                                                 Matrix_small_vector_column,
                                                 Matrix_intrusive_set_column>
                              > > > > > > >;
 
  struct Column_z2_settings {
    Column_z2_settings() : entryConstructor() {}
    Column_z2_settings([[maybe_unused]] Characteristic characteristic) : entryConstructor() {}
    Column_z2_settings([[maybe_unused]] const Column_z2_settings& toCopy) : entryConstructor() {}
    Column_z2_settings([[maybe_unused]] Column_z2_settings&& toMove) noexcept : entryConstructor() {}
    ~Column_z2_settings() = default;
    Column_z2_settings& operator=([[maybe_unused]] const Column_z2_settings& other) { return *this; }
    Column_z2_settings& operator=([[maybe_unused]] Column_z2_settings&& other) noexcept { return *this; }
 
    Entry_constructor entryConstructor;  // will be replaced by more specific allocators depending on the column type.
  };
 
  struct Column_zp_settings {
    Column_zp_settings() : operators(), entryConstructor() {}
    Column_zp_settings(Characteristic characteristic) : operators(characteristic), entryConstructor() {}
    Column_zp_settings(const Column_zp_settings& toCopy) : operators(toCopy.operators), entryConstructor() {}
    Column_zp_settings(Column_zp_settings&& toMove) noexcept
        : operators(std::move(toMove.operators)), entryConstructor() {}
    ~Column_zp_settings() = default;
    Column_zp_settings& operator=(const Column_zp_settings& other)
    {
      operators = other.operators;
      return *this;
    }
    Column_zp_settings& operator=(Column_zp_settings&& other) noexcept
    {
      if (this == &other) return *this;
 
      operators = std::move(other.operators);
      return *this;
    }
 
    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 = std::conditional_t<PersistenceMatrixOptions::is_z2, Column_z2_settings, Column_zp_settings>;

  static Field_operators const* get_operator_ptr(Column_settings const* colSettings)
  {
    if constexpr (PersistenceMatrixOptions::is_z2) {
      return nullptr;
    } else {
      if (colSettings == nullptr) return nullptr; // used for dummy columns
      return &(colSettings->operators);
    }
  }
  template <typename T>
  static Element get_coefficient_value(T v, [[maybe_unused]] Field_operators const* operators)
  {
    if constexpr (PersistenceMatrixOptions::is_z2) {
      return Field_operators::get_value(v);
    } else {
      return operators->get_value(v);
    }
  }
  static Element get_coefficient_value(bool v, [[maybe_unused]] Field_operators const* operators) { return v; }
 
  // 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 = std::conditional_t<PersistenceMatrixOptions::has_map_column_container,
                                              std::unordered_map<Index, Column>,
                                              std::vector<Column> >;
 
  static const bool hasFixedBarcode = Option_list::is_of_boundary_type && !PersistenceMatrixOptions::has_vine_update;
  using Barcode = std::conditional_t<
      hasFixedBarcode,
      std::vector<Bar>,
      std::conditional_t<PersistenceMatrixOptions::has_removable_columns, std::list<Bar>, std::vector<Bar> > >;
  using Bar_dictionary =
      std::conditional_t<hasFixedBarcode,
                         std::conditional_t<PersistenceMatrixOptions::can_retrieve_representative_cycles,
                                            std::vector<Index>,                   // RU
                                            std::unordered_map<Pos_index, Index>  // boundary
                                            >,
                         std::conditional_t<PersistenceMatrixOptions::has_removable_columns,
                                            std::unordered_map<Pos_index, typename Barcode::iterator>,
                                            std::vector<Index> > >;
 
  // default type for boundaries to permit list initialization directly in function parameters
  using Boundary = std::conditional_t<PersistenceMatrixOptions::is_z2,
                                      std::initializer_list<ID_index>,
                                      std::initializer_list<std::pair<ID_index, Element> > >;
 
  // 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 =
      std::conditional_t<PersistenceMatrixOptions::has_matrix_maximal_dimension_access || maxDimensionIsNeeded,
                         std::conditional_t<PersistenceMatrixOptions::has_removable_columns,
                                            Matrix_all_dimension_holder<Dimension>,
                                            Matrix_max_dimension_holder<Dimension> >,
                         Dummy_matrix_dimension_holder>;
 
  using Master_base_matrix = std::conditional_t<PersistenceMatrixOptions::has_column_compression,
                                                Base_matrix_with_column_compression<Matrix<PersistenceMatrixOptions> >,
                                                Base_matrix<Matrix<PersistenceMatrixOptions> > >;
  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 = std::conditional_t<PersistenceMatrixOptions::has_vine_update ||
                                                  PersistenceMatrixOptions::has_column_and_row_swaps,
                                              Base_swap<Matrix<PersistenceMatrixOptions>, Base>,
                                              Dummy_base_swap>;
  using Base_pairing_option =
      std::conditional_t<PersistenceMatrixOptions::has_column_pairings && !PersistenceMatrixOptions::has_vine_update &&
                             !PersistenceMatrixOptions::can_retrieve_representative_cycles,
                         Base_pairing<Matrix<PersistenceMatrixOptions> >,
                         Dummy_base_pairing>;
 
  using RU_pairing_option = std::conditional_t<PersistenceMatrixOptions::has_column_pairings,
                                               std::conditional_t<PersistenceMatrixOptions::has_vine_update,
                                                                  RU_barcode_swap<Matrix<PersistenceMatrixOptions> >,
                                                                  RU_pairing<Matrix<PersistenceMatrixOptions> > >,
                                               Dummy_ru_pairing>;
  using RU_vine_swap_option = std::conditional_t<PersistenceMatrixOptions::has_vine_update,
                                                 RU_vine_swap<Matrix<PersistenceMatrixOptions> >,
                                                 Dummy_ru_vine_swap>;
  using RU_representative_cycles_option =
      std::conditional_t<PersistenceMatrixOptions::can_retrieve_representative_cycles,
                         RU_representative_cycles<Matrix<PersistenceMatrixOptions> >,
                         Dummy_ru_representative_cycles>;
 
  using Chain_pairing_option =
      std::conditional_t<PersistenceMatrixOptions::has_column_pairings,
                         std::conditional_t<PersistenceMatrixOptions::has_vine_update,
                                            Chain_barcode_swap<Matrix<PersistenceMatrixOptions> >,
                                            Chain_pairing<Matrix<PersistenceMatrixOptions> > >,
                         Dummy_chain_pairing>;
  using Chain_vine_swap_option = std::conditional_t<PersistenceMatrixOptions::has_vine_update,
                                                    Chain_vine_swap<Matrix<PersistenceMatrixOptions> >,
                                                    Dummy_chain_vine_swap>;
  using Chain_representative_cycles_option =
      std::conditional_t<PersistenceMatrixOptions::can_retrieve_representative_cycles,
                         Chain_representative_cycles<Matrix<PersistenceMatrixOptions> >,
                         Dummy_chain_representative_cycles>;

  using Cycle = std::vector<Entry_representative>;

  template <class EntryRange>
  static Cycle build_cycle_from_range(const EntryRange& entries)
  {
    Cycle cycle;
    if constexpr (RangeTraits<EntryRange>::has_size) {
      cycle.reserve(entries.size());
    }
    for (const auto& c : entries) {
      if constexpr (PersistenceMatrixOptions::is_z2) {
        cycle.push_back(c.get_row_index());
      } else {
        cycle.push_back({c.get_row_index(), c.get_element()});
      }
    }
    return cycle;
  }
 
  // Return types to factorize the corresponding methods
 
  // The returned column is `const` if the matrix uses column compression
  using Returned_column =
      std::conditional_t<!isNonBasic && PersistenceMatrixOptions::has_column_compression, const Column, Column>;
  // The returned row is `const` if the matrix uses column compression
  using Returned_row =
      std::conditional_t<!isNonBasic && PersistenceMatrixOptions::has_column_compression, const Row, Row>;
  // 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 =
      std::conditional_t<PersistenceMatrixOptions::is_of_boundary_type || !isNonBasic ||
                             PersistenceMatrixOptions::column_indexation_type == Column_indexation_types::POSITION,
                         void,
                         std::vector<Entry_representative> >;

  Matrix();
  template <class Container = Boundary>
  Matrix(const std::vector<Container>& columns, Characteristic characteristic = 11);
  Matrix(unsigned int numberOfColumns, Characteristic characteristic = Field_operators::nullCharacteristic);
  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 = Field_operators::nullCharacteristic);
  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, class = std::enable_if_t<!std::is_arithmetic_v<Container> > >
  void insert_column(const Container& column);
  template <class Container, class = std::enable_if_t<!std::is_arithmetic_v<Container> > >
  void insert_column(const Container& column, Index columnIndex);
  void insert_column(ID_index idx, Element e = 1U);
  // 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>());

  template <class Boundary_range = Boundary>
  Insertion_return insert_maximal_cell(Index columnIndex,
                                       const Boundary_range& boundary,
                                       Dimension dim = Matrix::get_null_value<Dimension>());

  template <class Boundary_range = Boundary>
  Insertion_return insert_maximal_cell(Index columnIndex,
                                       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) &;
  Matrix& operator=(Matrix&& other) && noexcept;

  friend void swap(Matrix& matrix1, Matrix& matrix2) noexcept
  {
    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);

  void update_all_representative_cycles(Dimension dim = get_null_value<Dimension>());
  void update_representative_cycle(const Bar& bar);
  const std::vector<Cycle>& get_all_representative_cycles() const;
  const Cycle& get_representative_cycle(const Bar& bar) const;
 
 private:
  using Underlying_matrix = std::conditional_t<
      isNonBasic,
      std::conditional_t<
          PersistenceMatrixOptions::is_of_boundary_type,
          std::conditional_t<
              PersistenceMatrixOptions::has_vine_update || PersistenceMatrixOptions::can_retrieve_representative_cycles,
              std::conditional_t<
                  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>
              > >,
              std::conditional_t<
                  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>
              > > >,
          std::conditional_t<
              PersistenceMatrixOptions::column_indexation_type == Column_indexation_types::CONTAINER,
              Master_chain_matrix,
              std::conditional_t<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>
                                > > > >,
      Master_base_matrix>;
 
  // 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();
 
  Element _get_value(int coefficient) const;
};
 
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_);  // to avoid crashes at destruction, all columns have to be destroyed first
  delete colSettings_;
}
 
template <class PersistenceMatrixOptions>
inline void Matrix<PersistenceMatrixOptions>::set_characteristic(Characteristic characteristic)
{
  if constexpr (!PersistenceMatrixOptions::is_z2) {
    if (colSettings_->operators.get_characteristic() != Field_operators::nullCharacteristic) {
      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, class>
inline void Matrix<PersistenceMatrixOptions>::insert_column(const Container& column)
{
  if constexpr (!PersistenceMatrixOptions::is_z2) {
    GUDHI_CHECK(colSettings_->operators.get_characteristic() != Field_operators::nullCharacteristic,
                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, class>
inline void Matrix<PersistenceMatrixOptions>::insert_column(const Container& column, Index columnIndex)
{
  if constexpr (!PersistenceMatrixOptions::is_z2) {
    GUDHI_CHECK(colSettings_->operators.get_characteristic() != Field_operators::nullCharacteristic,
                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>
inline void Matrix<PersistenceMatrixOptions>::insert_column(ID_index idx, [[maybe_unused]] Element e)
{
  static_assert(
      !isNonBasic,
      "'insert_column' not available for the chosen options. The input has to be in the form of a cell boundary.");
 
  if constexpr (PersistenceMatrixOptions::is_z2) {
    matrix_.insert_column(idx);
  } else {
    GUDHI_CHECK(colSettings_->operators.get_characteristic() != Field_operators::nullCharacteristic,
                std::logic_error("Matrix::insert_column - Columns cannot be initialized if the coefficient field "
                                 "characteristic is not specified."));
    matrix_.insert_column(idx, e);
  }
}
 
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() != Field_operators::nullCharacteristic,
                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() != Field_operators::nullCharacteristic,
                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>
template <class Boundary_range>
inline typename Matrix<PersistenceMatrixOptions>::Insertion_return
Matrix<PersistenceMatrixOptions>::insert_maximal_cell(Index columnIndex, const Boundary_range& boundary, Dimension dim)
{
  static_assert(isNonBasic && PersistenceMatrixOptions::has_vine_update,
                "'insert_maximal_cell' is not available for the chosen options.");
 
  if constexpr (isNonBasic && !PersistenceMatrixOptions::is_of_boundary_type &&
                PersistenceMatrixOptions::column_indexation_type == Column_indexation_types::CONTAINER)
    return matrix_.insert_maximal_cell(columnIndex, boundary, dim);
  else
    matrix_.insert_maximal_cell(columnIndex, boundary, dim);
}
 
template <class PersistenceMatrixOptions>
template <class Boundary_range>
inline typename Matrix<PersistenceMatrixOptions>::Insertion_return
Matrix<PersistenceMatrixOptions>::insert_maximal_cell(Index columnIndex,
                                                      ID_index cellIndex,
                                                      const Boundary_range& boundary,
                                                      Dimension dim)
{
  static_assert(isNonBasic && PersistenceMatrixOptions::has_vine_update,
                "'insert_maximal_cell' is not available for the chosen options.");
 
  if constexpr (isNonBasic && !PersistenceMatrixOptions::is_of_boundary_type &&
                PersistenceMatrixOptions::column_indexation_type == Column_indexation_types::CONTAINER)
    return matrix_.insert_maximal_cell(columnIndex, cellIndex, boundary, dim);
  else
    matrix_.insert_maximal_cell(columnIndex, 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)
{
  matrix_.multiply_target_and_add_to(sourceColumnIndex, _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.");
 
  matrix_.multiply_target_and_add_to(sourceColumn, _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)
{
  matrix_.multiply_source_and_add_to(_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.");
 
  matrix_.multiply_source_and_add_to(_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 Matrix<PersistenceMatrixOptions>& Matrix<PersistenceMatrixOptions>::operator=(Matrix&& other) && noexcept
{
  if (this == &other) return *this;
 
  matrix_ = std::move(other.matrix_);
  colSettings_ = std::exchange(other.colSettings_, nullptr);
 
  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_all_representative_cycles(Dimension dim)
{
  static_assert(PersistenceMatrixOptions::can_retrieve_representative_cycles, "This method was not enabled.");
  matrix_.update_all_representative_cycles(dim);
}
 
template <class PersistenceMatrixOptions>
inline void Matrix<PersistenceMatrixOptions>::update_representative_cycle(const Bar& bar)
{
  static_assert(PersistenceMatrixOptions::can_retrieve_representative_cycles, "This method was not enabled.");
  matrix_.update_representative_cycle(bar);
}
 
template <class PersistenceMatrixOptions>
inline const std::vector<typename Matrix<PersistenceMatrixOptions>::Cycle>&
Matrix<PersistenceMatrixOptions>::get_all_representative_cycles() const
{
  static_assert(PersistenceMatrixOptions::can_retrieve_representative_cycles, "This method was not enabled.");
  return matrix_.get_all_representative_cycles();
}
 
template <class PersistenceMatrixOptions>
inline const typename Matrix<PersistenceMatrixOptions>::Cycle&
Matrix<PersistenceMatrixOptions>::get_representative_cycle(const Bar& bar) const
{
  static_assert(PersistenceMatrixOptions::can_retrieve_representative_cycles, "This method was not enabled.");
  return matrix_.get_representative_cycle(bar);
}
 
template <class PersistenceMatrixOptions>
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.");
}
 
template <class PersistenceMatrixOptions>
inline typename Matrix<PersistenceMatrixOptions>::Element Matrix<PersistenceMatrixOptions>::_get_value(
    int coefficient) const
{
  return get_coefficient_value(coefficient, get_operator_ptr(colSettings_));
}
 
}  // namespace persistence_matrix
}  // namespace Gudhi
 
#endif  // MASTER_MATRIX_H