Chain_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):
 *      - YYYY/MM Author: Description of the modification
 */

 
#ifndef PM_CHAIN_MATRIX_H
#define PM_CHAIN_MATRIX_H
 
#include <iostream>  //print() only
#include <set>
#include <map>
#include <stdexcept>
#include <vector>
#include <utility>    //std::swap, std::move & std::exchange
#include <algorithm>  //std::sort
 
#include <gudhi/Persistence_matrix/Position_to_index_overlay.h>  //friend
#include <gudhi/Persistence_matrix/Id_to_index_overlay.h>        //friend
#include <gudhi/Persistence_matrix/index_mapper.h>
 
namespace Gudhi {
namespace persistence_matrix {

template <class Master_matrix>
class Chain_matrix : public Master_matrix::Matrix_dimension_option,
                     public Master_matrix::Chain_pairing_option,
                     public Master_matrix::Chain_vine_swap_option,
                     public Master_matrix::Chain_representative_cycles_option,
                     public Master_matrix::Matrix_row_access_option,
                     protected std::conditional_t<
                         Master_matrix::Option_list::has_vine_update &&
                             (Master_matrix::Option_list::has_column_pairings ||
                              Master_matrix::Option_list::can_retrieve_representative_cycles),
                         Index_mapper<typename Master_matrix::template Dictionary<typename Master_matrix::Pos_index>>,
                         Dummy_index_mapper>
{
 private:
  using Dim_opt = typename Master_matrix::Matrix_dimension_option;
  using Pair_opt = typename Master_matrix::Chain_pairing_option;
  using Swap_opt = typename Master_matrix::Chain_vine_swap_option;
  using Rep_opt = typename Master_matrix::Chain_representative_cycles_option;
  using RA_opt = typename Master_matrix::Matrix_row_access_option;
 
  static constexpr bool hasPivotToPosMap_ =
      Master_matrix::Option_list::has_vine_update && (Master_matrix::Option_list::has_column_pairings ||
                                                      Master_matrix::Option_list::can_retrieve_representative_cycles);
 
  using Pivot_to_pos_mapper_opt =
      std::conditional_t<hasPivotToPosMap_,
                         Index_mapper<typename Master_matrix::template Dictionary<typename Master_matrix::Pos_index>>,
                         Dummy_index_mapper>;
 
 public:
  using Field_operators = typename Master_matrix::Field_operators;
  using Field_element = typename Master_matrix::Element;                     
  using Column = typename Master_matrix::Column;                             
  using Row = typename Master_matrix::Row;                                   
  using Entry = typename Master_matrix::Matrix_entry;                        
  using Entry_constructor = typename Master_matrix::Entry_constructor;       
  using Column_settings = typename Master_matrix::Column_settings;           
  using Boundary = typename Master_matrix::Boundary;                         
  using Entry_representative = typename Master_matrix::Entry_representative; 
  using Index = typename Master_matrix::Index;                               
  using ID_index = typename Master_matrix::ID_index;                         
  using Pos_index = typename Master_matrix::Pos_index;                       
  using Dimension = typename Master_matrix::Dimension;                       

  Chain_matrix(Column_settings* colSettings);
  template <class Boundary_range = Boundary>
  Chain_matrix(const std::vector<Boundary_range>& orderedBoundaries, Column_settings* colSettings);
  Chain_matrix(unsigned int numberOfColumns, Column_settings* colSettings);
  template <typename BirthComparatorFunction, typename DeathComparatorFunction>
  Chain_matrix(Column_settings* colSettings,
               const BirthComparatorFunction& birthComparator,
               const DeathComparatorFunction& deathComparator);
  template <typename BirthComparatorFunction, typename DeathComparatorFunction, class Boundary_range = Boundary>
  Chain_matrix(const std::vector<Boundary_range>& orderedBoundaries,
               Column_settings* colSettings,
               const BirthComparatorFunction& birthComparator,
               const DeathComparatorFunction& deathComparator);
  template <typename BirthComparatorFunction, typename DeathComparatorFunction>
  Chain_matrix(unsigned int numberOfColumns,
               Column_settings* colSettings,
               const BirthComparatorFunction& birthComparator,
               const DeathComparatorFunction& deathComparator);
  Chain_matrix(const Chain_matrix& matrixToCopy, Column_settings* colSettings = nullptr);
  Chain_matrix(Chain_matrix&& other) noexcept;
 
  ~Chain_matrix() = default;

  template <class Boundary_range = Boundary>
  std::vector<Entry_representative> insert_boundary(
      const Boundary_range& boundary,
      Dimension dim = Master_matrix::template get_null_value<Dimension>());
  template <class Boundary_range = Boundary>
  std::vector<Entry_representative> insert_boundary(
      ID_index cellID,
      const Boundary_range& boundary,
      Dimension dim = Master_matrix::template get_null_value<Dimension>());
  template <class Boundary_range = Boundary>
  std::vector<Entry_representative> insert_maximal_cell(
      Index columnIndex,
      const Boundary_range& boundary,
      Dimension dim = Master_matrix::template get_null_value<Dimension>());
  template <class Boundary_range = Boundary>
  std::vector<Entry_representative> insert_maximal_cell(
      Index columnIndex,
      ID_index cellID,
      const Boundary_range& boundary,
      Dimension dim = Master_matrix::template get_null_value<Dimension>());
  Column& get_column(Index columnIndex);
  const Column& get_column(Index columnIndex) const;
  void remove_maximal_cell(ID_index cellID);
  void remove_maximal_cell(ID_index cellID, const std::vector<ID_index>& columnsToSwap);
  void remove_last();

  Index get_number_of_columns() const;

  Dimension get_column_dimension(Index columnIndex) const;

  void add_to(Index sourceColumnIndex, Index targetColumnIndex);
  void multiply_target_and_add_to(Index sourceColumnIndex, const Field_element& coefficient, Index targetColumnIndex);
  void multiply_source_and_add_to(const Field_element& coefficient, Index sourceColumnIndex, Index targetColumnIndex);

  bool is_zero_entry(Index columnIndex, ID_index rowIndex) const;
  bool is_zero_column(Index columnIndex);

  Index get_column_with_pivot(ID_index cellID) const;
  ID_index get_pivot(Index columnIndex);

  void reset(Column_settings* colSettings)
  {
    if constexpr (Master_matrix::Option_list::has_matrix_maximal_dimension_access) Dim_opt::_reset();
    if constexpr (Master_matrix::Option_list::has_column_pairings) Pair_opt::_reset();
    if constexpr (Master_matrix::Option_list::can_retrieve_representative_cycles) Rep_opt::_reset();
    if constexpr (hasPivotToPosMap_) Pivot_to_pos_mapper_opt::map_.clear();
    matrix_.clear();
    pivotToColumnIndex_.clear();
    nextIndex_ = 0;
    nextPosition_ = 0;
    colSettings_ = colSettings;
  }

  Chain_matrix& operator=(const Chain_matrix& other);
  Chain_matrix& operator=(Chain_matrix&& other) noexcept;

  friend void swap(Chain_matrix& matrix1, Chain_matrix& matrix2) noexcept
  {
    swap(static_cast<Dim_opt&>(matrix1), static_cast<Dim_opt&>(matrix2));
    swap(static_cast<Pair_opt&>(matrix1), static_cast<Pair_opt&>(matrix2));
    swap(static_cast<Swap_opt&>(matrix1), static_cast<Swap_opt&>(matrix2));
    swap(static_cast<Rep_opt&>(matrix1), static_cast<Rep_opt&>(matrix2));
    swap(static_cast<Pivot_to_pos_mapper_opt&>(matrix1), static_cast<Pivot_to_pos_mapper_opt&>(matrix2));
    matrix1.matrix_.swap(matrix2.matrix_);
    matrix1.pivotToColumnIndex_.swap(matrix2.pivotToColumnIndex_);
    std::swap(matrix1.nextIndex_, matrix2.nextIndex_);
    std::swap(matrix1.nextPosition_, matrix2.nextPosition_);
    std::swap(matrix1.colSettings_, matrix2.colSettings_);
 
    if constexpr (Master_matrix::Option_list::has_row_access) {
      swap(static_cast<RA_opt&>(matrix1), static_cast<RA_opt&>(matrix2));
    }
  }
 
  void print() const;  // for debug
 
  friend class Id_to_index_overlay<Chain_matrix<Master_matrix>, Master_matrix>;
  friend class Position_to_index_overlay<Chain_matrix<Master_matrix>, Master_matrix>;
 
 private:
  using Column_container = typename Master_matrix::Column_container;
  using Dictionary = typename Master_matrix::template Dictionary<Index>;
  using Barcode = typename Master_matrix::Barcode;
  using Bar_dictionary = typename Master_matrix::Bar_dictionary;
  using Tmp_column = typename std::
      conditional<Master_matrix::Option_list::is_z2, std::set<ID_index>, std::map<ID_index, Field_element>>::type;
 
  friend Swap_opt;  // direct access to index mapper
  friend Pair_opt;  // direct access to index mapper
  friend Rep_opt;   // direct access to index mapper
 
  Column_container matrix_;       
  Dictionary pivotToColumnIndex_; 
  Index nextIndex_;               
  Pos_index nextPosition_;        
  Column_settings* colSettings_;  
 
  template <class Boundary_range>
  std::vector<Entry_representative> _reduce_boundary(ID_index cellID, const Boundary_range& boundary, Dimension dim);
  void _reduce_by_G(Tmp_column& column, std::vector<Entry_representative>& chainsInH, Index currentIndex);
  void _reduce_by_F(Tmp_column& column, std::vector<Entry_representative>& chainsInF, Index currentIndex);
  void _build_from_H(ID_index cellID, Tmp_column& column, std::vector<Entry_representative>& chainsInH);
  void _update_largest_death_in_F(const std::vector<Entry_representative>& chainsInF);
  void _insert_chain(const Tmp_column& column, Dimension dimension);
  void _insert_chain(const Tmp_column& column, Dimension dimension, Index pair);
  void _add_to(const Column& column, Tmp_column& set, unsigned int coef);
  template <typename F>
  void _add_to(Column& target, F&& addition);
  void _remove_last(Index lastIndex);
  void _update_barcode(Pos_index birth);
  void _add_bar(Dimension dim);
  template <class Container>
  void _container_insert(const Container& column, Index pos, Dimension dim);
  template <class ColumnIterator>
  void _container_insert(const ColumnIterator& rep);
 
  static void _insert_in(ID_index cellID, Tmp_column& column);
};
 
template <class Master_matrix>
inline Chain_matrix<Master_matrix>::Chain_matrix(Column_settings* colSettings)
    : Dim_opt(Master_matrix::template get_null_value<Dimension>()),
      Pair_opt(),
      Swap_opt(),
      Rep_opt(),
      RA_opt(),
      Pivot_to_pos_mapper_opt(),
      nextIndex_(0),
      nextPosition_(0),
      colSettings_(colSettings)
{}
 
template <class Master_matrix>
template <class Boundary_range>
inline Chain_matrix<Master_matrix>::Chain_matrix(const std::vector<Boundary_range>& orderedBoundaries,
                                                 Column_settings* colSettings)
    : Dim_opt(Master_matrix::template get_null_value<Dimension>()),
      Pair_opt(),
      Swap_opt(),
      Rep_opt(),
      RA_opt(orderedBoundaries.size()),
      Pivot_to_pos_mapper_opt(),
      nextIndex_(0),
      nextPosition_(0),
      colSettings_(colSettings)
{
  matrix_.reserve(orderedBoundaries.size());
  if constexpr (Master_matrix::Option_list::has_map_column_container) {
    pivotToColumnIndex_.reserve(orderedBoundaries.size());
  } else {
    pivotToColumnIndex_.resize(orderedBoundaries.size(), Master_matrix::template get_null_value<Index>());
  }
 
  for (const Boundary_range& b : orderedBoundaries) {
    insert_boundary(b);
  }
}
 
template <class Master_matrix>
inline Chain_matrix<Master_matrix>::Chain_matrix(unsigned int numberOfColumns, Column_settings* colSettings)
    : Dim_opt(Master_matrix::template get_null_value<Dimension>()),
      Pair_opt(),
      Swap_opt(),
      Rep_opt(),
      RA_opt(numberOfColumns),
      Pivot_to_pos_mapper_opt(),
      nextIndex_(0),
      nextPosition_(0),
      colSettings_(colSettings)
{
  matrix_.reserve(numberOfColumns);
  if constexpr (Master_matrix::Option_list::has_map_column_container) {
    pivotToColumnIndex_.reserve(numberOfColumns);
  } else {
    pivotToColumnIndex_.resize(numberOfColumns, Master_matrix::template get_null_value<Index>());
  }
}
 
template <class Master_matrix>
template <typename BirthComparatorFunction, typename DeathComparatorFunction>
inline Chain_matrix<Master_matrix>::Chain_matrix(Column_settings* colSettings,
                                                 const BirthComparatorFunction& birthComparator,
                                                 const DeathComparatorFunction& deathComparator)
    : Dim_opt(Master_matrix::template get_null_value<Dimension>()),
      Pair_opt(),
      Swap_opt(birthComparator, deathComparator),
      Rep_opt(),
      RA_opt(),
      Pivot_to_pos_mapper_opt(),
      nextIndex_(0),
      nextPosition_(0),
      colSettings_(colSettings)
{}
 
template <class Master_matrix>
template <typename BirthComparatorFunction, typename DeathComparatorFunction, class Boundary_range>
inline Chain_matrix<Master_matrix>::Chain_matrix(const std::vector<Boundary_range>& orderedBoundaries,
                                                 Column_settings* colSettings,
                                                 const BirthComparatorFunction& birthComparator,
                                                 const DeathComparatorFunction& deathComparator)
    : Dim_opt(Master_matrix::template get_null_value<Dimension>()),
      Pair_opt(),
      Swap_opt(birthComparator, deathComparator),
      Rep_opt(),
      RA_opt(orderedBoundaries.size()),
      Pivot_to_pos_mapper_opt(),
      nextIndex_(0),
      nextPosition_(0),
      colSettings_(colSettings)
{
  matrix_.reserve(orderedBoundaries.size());
  if constexpr (Master_matrix::Option_list::has_map_column_container) {
    pivotToColumnIndex_.reserve(orderedBoundaries.size());
  } else {
    pivotToColumnIndex_.resize(orderedBoundaries.size(), Master_matrix::template get_null_value<Index>());
  }
  for (const Boundary_range& b : orderedBoundaries) {
    insert_boundary(b);
  }
}
 
template <class Master_matrix>
template <typename BirthComparatorFunction, typename DeathComparatorFunction>
inline Chain_matrix<Master_matrix>::Chain_matrix(unsigned int numberOfColumns,
                                                 Column_settings* colSettings,
                                                 const BirthComparatorFunction& birthComparator,
                                                 const DeathComparatorFunction& deathComparator)
    : Dim_opt(Master_matrix::template get_null_value<Dimension>()),
      Pair_opt(),
      Swap_opt(birthComparator, deathComparator),
      Rep_opt(),
      RA_opt(numberOfColumns),
      Pivot_to_pos_mapper_opt(),
      nextIndex_(0),
      nextPosition_(0),
      colSettings_(colSettings)
{
  matrix_.reserve(numberOfColumns);
  if constexpr (Master_matrix::Option_list::has_map_column_container) {
    pivotToColumnIndex_.reserve(numberOfColumns);
  } else {
    pivotToColumnIndex_.resize(numberOfColumns, Master_matrix::template get_null_value<Index>());
  }
}
 
template <class Master_matrix>
inline Chain_matrix<Master_matrix>::Chain_matrix(const Chain_matrix& matrixToCopy, Column_settings* colSettings)
    : Dim_opt(static_cast<const Dim_opt&>(matrixToCopy)),
      Pair_opt(static_cast<const Pair_opt&>(matrixToCopy)),
      Swap_opt(static_cast<const Swap_opt&>(matrixToCopy)),
      Rep_opt(static_cast<const Rep_opt&>(matrixToCopy)),
      RA_opt(static_cast<const RA_opt&>(matrixToCopy)),
      Pivot_to_pos_mapper_opt(static_cast<const Pivot_to_pos_mapper_opt&>(matrixToCopy)),
      pivotToColumnIndex_(matrixToCopy.pivotToColumnIndex_),
      nextIndex_(matrixToCopy.nextIndex_),
      nextPosition_(matrixToCopy.nextPosition_),
      colSettings_(colSettings == nullptr ? matrixToCopy.colSettings_ : colSettings)
{
  matrix_.reserve(matrixToCopy.matrix_.size());
  for (const auto& cont : matrixToCopy.matrix_) {
    _container_insert(cont);
  }
}
 
template <class Master_matrix>
inline Chain_matrix<Master_matrix>::Chain_matrix(Chain_matrix&& other) noexcept
    : Dim_opt(std::move(static_cast<Dim_opt&>(other))),
      Pair_opt(std::move(static_cast<Pair_opt&>(other))),
      Swap_opt(std::move(static_cast<Swap_opt&>(other))),
      Rep_opt(std::move(static_cast<Rep_opt&>(other))),
      RA_opt(std::move(static_cast<RA_opt&>(other))),
      Pivot_to_pos_mapper_opt(std::move(static_cast<Pivot_to_pos_mapper_opt&>(other))),
      matrix_(std::move(other.matrix_)),
      pivotToColumnIndex_(std::move(other.pivotToColumnIndex_)),
      nextIndex_(std::exchange(other.nextIndex_, 0)),
      nextPosition_(std::exchange(other.nextPosition_, 0)),
      colSettings_(std::exchange(other.colSettings_, nullptr))
{}
 
template <class Master_matrix>
template <class Boundary_range>
inline std::vector<typename Master_matrix::Entry_representative> Chain_matrix<Master_matrix>::insert_boundary(
    const Boundary_range& boundary,
    Dimension dim)
{
  return insert_boundary(nextIndex_, boundary, dim);
}
 
template <class Master_matrix>
template <class Boundary_range>
inline std::vector<typename Master_matrix::Entry_representative>
Chain_matrix<Master_matrix>::insert_boundary(ID_index cellID, const Boundary_range& boundary, Dimension dim)
{
  if constexpr (!Master_matrix::Option_list::has_map_column_container) {
    if (pivotToColumnIndex_.size() <= cellID) {
      pivotToColumnIndex_.resize((cellID * 2) + 1, Master_matrix::template get_null_value<Index>());
    }
  }
 
  if constexpr (hasPivotToPosMap_) {
    if constexpr (Master_matrix::Option_list::has_map_column_container) {
      Pivot_to_pos_mapper_opt::map_.try_emplace(cellID, nextPosition_);
    } else {
      if (Pivot_to_pos_mapper_opt::map_.size() <= cellID)
        Pivot_to_pos_mapper_opt::map_.resize(pivotToColumnIndex_.size(),
                                             Master_matrix::template get_null_value<Pos_index>());
      Pivot_to_pos_mapper_opt::map_[cellID] = nextPosition_;
    }
  }
 
  if constexpr (Master_matrix::Option_list::has_matrix_maximal_dimension_access) {
    Dim_opt::_update_up(dim == Master_matrix::template get_null_value<Dimension>()
                            ? (boundary.size() == 0 ? 0 : boundary.size() - 1)
                            : dim);
  }
 
  return _reduce_boundary(cellID, boundary, dim);
}
 
template <class Master_matrix>
template <class Boundary_range>
inline std::vector<typename Master_matrix::Entry_representative>
Chain_matrix<Master_matrix>::insert_maximal_cell(Index columnIndex, const Boundary_range& boundary, Dimension dim)
{
  return insert_maximal_cell(columnIndex, nextIndex_, boundary, dim);
}
 
template <class Master_matrix>
template <class Boundary_range>
inline std::vector<typename Master_matrix::Entry_representative> Chain_matrix<Master_matrix>::insert_maximal_cell(
    Index columnIndex,
    ID_index cellID,
    const Boundary_range& boundary,
    Dimension dim)
{
  static_assert(
      Master_matrix::Option_list::has_vine_update && (Master_matrix::Option_list::has_column_pairings ||
                                                      Master_matrix::Option_list::can_retrieve_representative_cycles),
      "'insert_maximal_cell' is not implemented for the chosen options.");
 
  GUDHI_CHECK(columnIndex >= 0, std::invalid_argument("Indices have to be positive."));
 
  auto chainsInF = insert_boundary(cellID, boundary, dim);
 
  const auto& pivotToPosition = Pivot_to_pos_mapper_opt::map_;
  Pos_index startPos = nextPosition_ - 1;
  Index startIndex = pivotToColumnIndex_[cellID];
 
  if (startPos > columnIndex) {
    std::vector<Index> colToSwap;
    colToSwap.reserve(pivotToPosition.size());
 
    // first loop stores all column indices whose filtration position should be after the new column
    // then orders the matrix indices by filtration position as they do not have to correspond
    // the underlying column container is not the same depending on `has_map_column_container`: one is a map, the
    // other a vector. To optimize the usage of both containers, the `if constexpr` differentiation is made
    if constexpr (Master_matrix::Option_list::has_map_column_container) {
      for (auto& p : pivotToPosition) {
        if (p.second >= columnIndex && p.second < startPos) colToSwap.push_back(pivotToColumnIndex_.at(p.first));
      }
      std::sort(colToSwap.begin(), colToSwap.end(), [&](Index c1, Index c2) {
        return pivotToPosition.at(get_pivot(c1)) > pivotToPosition.at(get_pivot(c2));
      });
    } else {
      for (Index i = 0; i < pivotToPosition.size(); ++i) {
        if (pivotToPosition[i] >= columnIndex && pivotToPosition[i] < startPos)
          colToSwap.push_back(pivotToColumnIndex_[i]);
      }
      std::sort(colToSwap.begin(), colToSwap.end(), [&](Index c1, Index c2) {
        return pivotToPosition[get_pivot(c1)] > pivotToPosition[get_pivot(c2)];
      });
    }
 
    for (Index i : colToSwap) {
      startIndex = Swap_opt::vine_swap(i, startIndex) == startIndex ? i : startIndex;
    }
  }
 
  return chainsInF;
}
 
template <class Master_matrix>
inline typename Chain_matrix<Master_matrix>::Column& Chain_matrix<Master_matrix>::get_column(Index columnIndex)
{
  if constexpr (Master_matrix::Option_list::has_map_column_container) {
    return matrix_.at(columnIndex);
  } else {
    return matrix_[columnIndex];
  }
}
 
template <class Master_matrix>
inline const typename Chain_matrix<Master_matrix>::Column& Chain_matrix<Master_matrix>::get_column(
    Index columnIndex) const
{
  if constexpr (Master_matrix::Option_list::has_map_column_container) {
    return matrix_.at(columnIndex);
  } else {
    return matrix_[columnIndex];
  }
}
 
template <class Master_matrix>
inline void Chain_matrix<Master_matrix>::remove_maximal_cell(ID_index cellID)
{
  static_assert(
      Master_matrix::Option_list::has_removable_columns && Master_matrix::Option_list::has_map_column_container,
      "'remove_maximal_cell' is not implemented for the chosen options.");
  static_assert(
      Master_matrix::Option_list::has_vine_update && (Master_matrix::Option_list::has_column_pairings ||
                                                      Master_matrix::Option_list::can_retrieve_representative_cycles),
      "'remove_maximal_cell' is not implemented for the chosen options.");
 
  // TODO: find simple test to verify that col at columnIndex is maximal even without row access.
 
  const auto& pivotToPosition = Pivot_to_pos_mapper_opt::map_;
  auto it = pivotToPosition.find(cellID);
  if (it == pivotToPosition.end()) return;  // cell does not exists. TODO: put an assert instead?
  Pos_index startPos = it->second;
  Index startIndex = pivotToColumnIndex_.at(cellID);
 
  if (startPos != nextPosition_ - 1) {
    std::vector<Index> colToSwap;
    colToSwap.reserve(matrix_.size());
 
    for (auto& p : pivotToPosition) {
      if (p.second > startPos) colToSwap.push_back(pivotToColumnIndex_.at(p.first));
    }
    std::sort(colToSwap.begin(), colToSwap.end(), [&](Index c1, Index c2) {
      return pivotToPosition.at(get_pivot(c1)) < pivotToPosition.at(get_pivot(c2));
    });
 
    for (Index i : colToSwap) {
      startIndex = Swap_opt::vine_swap(startIndex, i);
    }
  }
 
  _remove_last(startIndex);
}
 
template <class Master_matrix>
inline void Chain_matrix<Master_matrix>::remove_maximal_cell(ID_index cellID,
                                                             const std::vector<ID_index>& columnsToSwap)
{
  static_assert(Master_matrix::Option_list::has_removable_columns,
                "'remove_maximal_cell' is not implemented for the chosen options.");
  static_assert(Master_matrix::Option_list::has_map_column_container && Master_matrix::Option_list::has_vine_update,
                "'remove_maximal_cell' is not implemented for the chosen options.");
 
  // TODO: find simple test to verify that col at columnIndex is maximal even without row access.
 
  Index startIndex = pivotToColumnIndex_.at(cellID);
 
  for (ID_index i : columnsToSwap) {
    startIndex = Swap_opt::vine_swap(startIndex, pivotToColumnIndex_.at(i));
  }
 
  _remove_last(startIndex);
}
 
template <class Master_matrix>
inline void Chain_matrix<Master_matrix>::remove_last()
{
  static_assert(Master_matrix::Option_list::has_removable_columns,
                "'remove_last' is not implemented for the chosen options.");
  static_assert(Master_matrix::Option_list::has_map_column_container || !Master_matrix::Option_list::has_vine_update,
                "'remove_last' is not implemented for the chosen options.");
 
  if (nextIndex_ == 0 || matrix_.empty()) return;  // empty matrix
 
  if constexpr (Master_matrix::Option_list::has_vine_update) {
    // careful: linear because of the search of the last index. It is better to keep track of the @ref IDIdx index
    // of the last column while performing swaps (or the @ref MatIdx with the return values of `vine_swap` + get_pivot)
    // and then call `remove_maximal_cell` with it and an empty `columnsToSwap`.
 
    ID_index pivot = 0;
    Index colIndex = 0;
    for (auto& p : pivotToColumnIndex_) {
      if (p.first > pivot) {  // pivots have to be strictly increasing in order of filtration
        pivot = p.first;
        colIndex = p.second;
      }
    }
    _remove_last(colIndex);
  } else {
    _remove_last(nextIndex_ - 1);
  }
}
 
template <class Master_matrix>
inline typename Chain_matrix<Master_matrix>::Index Chain_matrix<Master_matrix>::get_number_of_columns() const
{
  if constexpr (Master_matrix::Option_list::has_map_column_container) {
    return matrix_.size();
  } else {
    return nextIndex_;  // matrix could have been resized much bigger while insert
  }
}
 
template <class Master_matrix>
inline typename Chain_matrix<Master_matrix>::Dimension Chain_matrix<Master_matrix>::get_column_dimension(
    Index columnIndex) const
{
  return get_column(columnIndex).get_dimension();
}
 
template <class Master_matrix>
inline void Chain_matrix<Master_matrix>::add_to(Index sourceColumnIndex, Index targetColumnIndex)
{
  auto& col = get_column(targetColumnIndex);
  _add_to(col, [&]() { col += get_column(sourceColumnIndex); });
}
 
template <class Master_matrix>
inline void Chain_matrix<Master_matrix>::multiply_target_and_add_to(Index sourceColumnIndex,
                                                                    const Field_element& coefficient,
                                                                    Index targetColumnIndex)
{
  auto& col = get_column(targetColumnIndex);
  _add_to(col, [&]() { col.multiply_target_and_add(coefficient, get_column(sourceColumnIndex)); });
}
 
template <class Master_matrix>
inline void Chain_matrix<Master_matrix>::multiply_source_and_add_to(const Field_element& coefficient,
                                                                    Index sourceColumnIndex,
                                                                    Index targetColumnIndex)
{
  auto& col = get_column(targetColumnIndex);
  _add_to(col, [&]() { col.multiply_source_and_add(get_column(sourceColumnIndex), coefficient); });
}
 
template <class Master_matrix>
inline bool Chain_matrix<Master_matrix>::is_zero_entry(Index columnIndex, ID_index rowIndex) const
{
  return !get_column(columnIndex).is_non_zero(rowIndex);
}
 
template <class Master_matrix>
inline bool Chain_matrix<Master_matrix>::is_zero_column(Index columnIndex)
{
  return get_column(columnIndex).is_empty();
}
 
template <class Master_matrix>
inline typename Chain_matrix<Master_matrix>::Index Chain_matrix<Master_matrix>::get_column_with_pivot(
    ID_index cellID) const
{
  if constexpr (Master_matrix::Option_list::has_map_column_container) {
    return pivotToColumnIndex_.at(cellID);
  } else {
    return pivotToColumnIndex_[cellID];
  }
}
 
template <class Master_matrix>
inline typename Chain_matrix<Master_matrix>::ID_index Chain_matrix<Master_matrix>::get_pivot(Index columnIndex)
{
  return get_column(columnIndex).get_pivot();
}
 
template <class Master_matrix>
inline Chain_matrix<Master_matrix>& Chain_matrix<Master_matrix>::operator=(const Chain_matrix& other)
{
  if (this == &other) return *this;
 
  Dim_opt::operator=(other);
  Swap_opt::operator=(other);
  Pair_opt::operator=(other);
  Rep_opt::operator=(other);
  Pivot_to_pos_mapper_opt::operator=(other);
  matrix_.clear();
  pivotToColumnIndex_ = other.pivotToColumnIndex_;
  nextIndex_ = other.nextIndex_;
  nextPosition_ = other.nextPosition_;
  colSettings_ = other.colSettings_;
 
  matrix_.reserve(other.matrix_.size());
  for (const auto& cont : other.matrix_) {
    _container_insert(cont);
  }
 
  return *this;
}
 
template <class Master_matrix>
inline Chain_matrix<Master_matrix>& Chain_matrix<Master_matrix>::operator=(Chain_matrix&& other) noexcept
{
  if (this == &other) return *this;
 
  Dim_opt::operator=(std::move(other));
  Pair_opt::operator=(std::move(other));
  Swap_opt::operator=(std::move(other));
  Rep_opt::operator=(std::move(other));
  RA_opt::operator=(std::move(other));
  Pivot_to_pos_mapper_opt::operator=(std::move(other));
  matrix_ = std::move(other.matrix_);
  pivotToColumnIndex_ = std::move(other.pivotToColumnIndex_);
  nextIndex_ = std::exchange(other.nextIndex_, 0);
  nextPosition_ = std::exchange(other.nextPosition_, 0);
  colSettings_ = std::exchange(other.colSettings_, nullptr);
 
  return *this;
}
 
template <class Master_matrix>
inline void Chain_matrix<Master_matrix>::print() const
{
  std::cout << "Column Matrix:\n";
  if constexpr (!Master_matrix::Option_list::has_map_column_container) {
    for (ID_index i = 0; i < pivotToColumnIndex_.size(); ++i) {
      Index pos = pivotToColumnIndex_[i];
      if (pos != Master_matrix::template get_null_value<Index>()) {
        const Column& col = matrix_[pos];
        for (const auto& entry : col) {
          std::cout << entry.get_row_index() << " ";
        }
        std::cout << "(" << i << ", " << pos << ")\n";
      }
    }
    if constexpr (Master_matrix::Option_list::has_row_access) {
      std::cout << "\n";
      std::cout << "Row Matrix:\n";
      for (ID_index i = 0; i < pivotToColumnIndex_.size(); ++i) {
        Index pos = pivotToColumnIndex_[i];
        if (pos != Master_matrix::template get_null_value<Index>()) {
          const Row& row = RA_opt::get_row(pos);
          for (const auto& entry : row) {
            std::cout << entry.get_column_index() << " ";
          }
          std::cout << "(" << i << ", " << pos << ")\n";
        }
      }
    }
  } else {
    for (const auto& p : pivotToColumnIndex_) {
      const Column& col = matrix_.at(p.second);
      for (const auto& entry : col) {
        std::cout << entry.get_row_index() << " ";
      }
      std::cout << "(" << p.first << ", " << p.second << ")\n";
    }
    if constexpr (Master_matrix::Option_list::has_row_access) {
      std::cout << "\n";
      std::cout << "Row Matrix:\n";
      for (const auto& p : pivotToColumnIndex_) {
        const Row& row = RA_opt::get_row(p.first);
        for (const auto& entry : row) {
          std::cout << entry.get_column_index() << " ";
        }
        std::cout << "(" << p.first << ", " << p.second << ")\n";
      }
    }
  }
  std::cout << "\n";
}
 
template <class Master_matrix>
template <class Boundary_range>
inline std::vector<typename Master_matrix::Entry_representative>
Chain_matrix<Master_matrix>::_reduce_boundary(ID_index cellID, const Boundary_range& boundary, Dimension dim)
{
  Tmp_column column(boundary.begin(), boundary.end());
  if (dim == Master_matrix::template get_null_value<Dimension>())
    dim = boundary.begin() == boundary.end() ? 0 : boundary.size() - 1;
  std::vector<Entry_representative> chainsInH;  // for corresponding indices in H (paired columns)
  std::vector<Entry_representative> chainsInF;  // for corresponding indices in F (unpaired, essential columns)
 
  auto get_last = [&column]() {
    if constexpr (Master_matrix::Option_list::is_z2)
      return *(column.rbegin());
    else
      return column.rbegin()->first;
  };
 
  if (boundary.begin() == boundary.end()) {
    _insert_in(cellID, column);
    _insert_chain(column, dim);
    return chainsInF;
  }
 
  Index currentIndex = get_column_with_pivot(get_last());
 
  while (get_column(currentIndex).is_paired()) {
    _reduce_by_G(column, chainsInH, currentIndex);
 
    if (column.empty()) {
      // produce the sum of all col_h in chains_in_H
      _build_from_H(cellID, column, chainsInH);
      // create a new cycle (in F) sigma - \sum col_h
      _insert_chain(column, dim);
      return chainsInF;
    }
 
    currentIndex = get_column_with_pivot(get_last());
  }
 
  while (!column.empty()) {
    currentIndex = get_column_with_pivot(get_last());
 
    if (!get_column(currentIndex).is_paired()) {
      // only fills currentEssentialCycleIndices if Z2 coefficients, so chainsInF remains empty
      _reduce_by_F(column, chainsInF, currentIndex);
    } else {
      _reduce_by_G(column, chainsInH, currentIndex);
    }
  }
 
  _update_largest_death_in_F(chainsInF);
 
  // Compute the new column cellID + \sum col_h, for col_h in chains_in_H
  _build_from_H(cellID, column, chainsInH);
 
  // Create and insert (\sum col_h) + sigma (in H, paired with chain_fp) in matrix_
  _insert_chain(column, dim, Master_matrix::get_row_index(chainsInF[0]));
 
  return chainsInF;
}
 
template <class Master_matrix>
inline void Chain_matrix<Master_matrix>::_reduce_by_G(Tmp_column& column,
                                                      std::vector<Entry_representative>& chainsInH,
                                                      Index currentIndex)
{
  Column& col = get_column(currentIndex);
  if constexpr (Master_matrix::Option_list::is_z2) {
    _add_to(col, column, 1U);                           // Reduce with the column col_g
    chainsInH.push_back(col.get_paired_chain_index());  // keep the col_h with which col_g is paired
  } else {
    Field_element coef = col.get_pivot_value();
    auto& operators = colSettings_->operators;
    coef = operators.get_inverse(coef);
    operators.multiply_inplace(coef, operators.get_characteristic() - column.rbegin()->second);
 
    _add_to(col, column, coef);                                  // Reduce with the column col_g
    chainsInH.emplace_back(col.get_paired_chain_index(), coef);  // keep the col_h with which col_g is paired
  }
}
 
template <class Master_matrix>
inline void Chain_matrix<Master_matrix>::_reduce_by_F(Tmp_column& column,
                                                      std::vector<Entry_representative>& chainsInF,
                                                      Index currentIndex)
{
  Column& col = get_column(currentIndex);
  if constexpr (Master_matrix::Option_list::is_z2) {
    _add_to(col, column, 1U);  // Reduce with the column col_g
    chainsInF.push_back(currentIndex);
  } else {
    Field_element coef = col.get_pivot_value();
    auto& operators = colSettings_->operators;
    coef = operators.get_inverse(coef);
    operators.multiply_inplace(coef, operators.get_characteristic() - column.rbegin()->second);
 
    _add_to(col, column, coef);  // Reduce with the column col_g
    chainsInF.emplace_back(currentIndex, operators.get_characteristic() - coef);
  }
}
 
template <class Master_matrix>
inline void Chain_matrix<Master_matrix>::_build_from_H(ID_index cellID,
                                                       Tmp_column& column,
                                                       std::vector<Entry_representative>& chainsInH)
{
  _insert_in(cellID, column);
 
  for (const auto& idx_h : chainsInH) {
    _add_to(get_column(Master_matrix::get_row_index(idx_h)), column, Master_matrix::get_element(idx_h));
  }
}
 
template <class Master_matrix>
inline void Chain_matrix<Master_matrix>::_update_largest_death_in_F(const std::vector<Entry_representative>& chainsInF)
{
  Index toUpdate = Master_matrix::get_row_index(chainsInF[0]);
  if constexpr (Master_matrix::Option_list::is_z2) {
    for (auto other_col_it = chainsInF.begin() + 1; other_col_it != chainsInF.end(); ++other_col_it) {
      add_to(*other_col_it, toUpdate);
    }
  } else {
    get_column(toUpdate) *= Master_matrix::get_element(chainsInF[0]);
    for (auto other_col_it = chainsInF.begin() + 1; other_col_it != chainsInF.end(); ++other_col_it) {
      multiply_source_and_add_to(
          Master_matrix::get_element(*other_col_it), Master_matrix::get_row_index(*other_col_it), toUpdate);
    }
  }
}
 
template <class Master_matrix>
inline void Chain_matrix<Master_matrix>::_insert_chain(const Tmp_column& column, Dimension dimension)
{
  _container_insert(column, nextIndex_, dimension);
  _add_bar(dimension);
 
  ++nextIndex_;
}
 
template <class Master_matrix>
inline void Chain_matrix<Master_matrix>::_insert_chain(const Tmp_column& column, Dimension dimension, Index pair)
{
  // true when no vine updates and if nextIndex_ is updated in remove_last for special case of no vines
  // because then @ref PosIdx == @ref MatIdx
  Pos_index pairPos = pair;
 
  _container_insert(column, nextIndex_, dimension);
 
  get_column(nextIndex_).assign_paired_chain(pair);
  auto& pairCol = get_column(pair);
  pairCol.assign_paired_chain(nextIndex_);
 
  if constexpr (Master_matrix::Option_list::has_column_pairings && Master_matrix::Option_list::has_vine_update) {
    pairPos = Pivot_to_pos_mapper_opt::map_[pairCol.get_pivot()];
  }
 
  _update_barcode(pairPos);
 
  ++nextIndex_;
}
 
template <class Master_matrix>
inline void Chain_matrix<Master_matrix>::_add_to(const Column& column,
                                                 Tmp_column& set,
                                                 [[maybe_unused]] unsigned int coef)
{
  if constexpr (Master_matrix::Option_list::is_z2) {
    std::pair<typename std::set<Index>::iterator, bool> res_insert;
    for (const Entry& entry : column) {
      res_insert = set.insert(entry.get_row_index());
      if (!res_insert.second) {
        set.erase(res_insert.first);
      }
    }
  } else {
    auto& operators = colSettings_->operators;
    for (const Entry& entry : column) {
      auto res = set.emplace(entry.get_row_index(), entry.get_element());
      if (res.second) {
        operators.multiply_inplace(res.first->second, coef);
      } else {
        operators.multiply_and_add_inplace_back(entry.get_element(), coef, res.first->second);
        if (res.first->second == Field_operators::get_additive_identity()) {
          set.erase(res.first);
        }
      }
    }
  }
}
 
template <class Master_matrix>
template <typename F>
inline void Chain_matrix<Master_matrix>::_add_to(Column& target, F&& addition)
{
  auto pivot = target.get_pivot();
  std::forward<F>(addition)();
 
  if (pivot != target.get_pivot()) {
    if constexpr (Master_matrix::Option_list::has_map_column_container) {
      std::swap(pivotToColumnIndex_.at(pivot), pivotToColumnIndex_.at(target.get_pivot()));
    } else {
      std::swap(pivotToColumnIndex_[pivot], pivotToColumnIndex_[target.get_pivot()]);
    }
  }
}
 
template <class Master_matrix>
inline void Chain_matrix<Master_matrix>::_remove_last(Index lastIndex)
{
  static_assert(Master_matrix::Option_list::has_removable_columns,
                "'_remove_last' is not implemented for the chosen options.");
  static_assert(Master_matrix::Option_list::has_map_column_container || !Master_matrix::Option_list::has_vine_update,
                "'_remove_last' is not implemented for the chosen options.");
 
  ID_index pivot;
 
  if constexpr (Master_matrix::Option_list::has_map_column_container) {
    auto itToErase = matrix_.find(lastIndex);
    Column& colToErase = itToErase->second;
    pivot = colToErase.get_pivot();
 
    if constexpr (Master_matrix::Option_list::has_matrix_maximal_dimension_access) {
      Dim_opt::_update_down(colToErase.get_dimension());
    }
 
    if (colToErase.is_paired()) matrix_.at(colToErase.get_paired_chain_index()).unassign_paired_chain();
    pivotToColumnIndex_.erase(pivot);
    matrix_.erase(itToErase);
  } else {
    GUDHI_CHECK(lastIndex == nextIndex_ - 1 && nextIndex_ == matrix_.size(),
                std::logic_error("Chain_matrix::_remove_last - Indexation problem."));
 
    Column& colToErase = matrix_[lastIndex];
    pivot = colToErase.get_pivot();
 
    if constexpr (Master_matrix::Option_list::has_matrix_maximal_dimension_access) {
      Dim_opt::_update_down(colToErase.get_dimension());
    }
 
    if (colToErase.is_paired()) matrix_.at(colToErase.get_paired_chain_index()).unassign_paired_chain();
    pivotToColumnIndex_[pivot] = Master_matrix::template get_null_value<Index>();
    matrix_.pop_back();
    // TODO: resize matrix_ when a lot is removed? Could be not the best strategy if user inserts a lot back afterwards.
  }
 
  if constexpr (!Master_matrix::Option_list::has_vine_update) {
    --nextIndex_;  // should not be updated when there are vine updates, as possibly lastIndex != nextIndex - 1
  }
 
  --nextPosition_;
  if constexpr (Master_matrix::Option_list::has_column_pairings) {
    Pair_opt::_erase_bar(nextPosition_);
  }
  if constexpr (hasPivotToPosMap_) {
    Pivot_to_pos_mapper_opt::map_.erase(pivot);
  }
 
  if constexpr (Master_matrix::Option_list::has_row_access) {
    GUDHI_CHECK(
        RA_opt::get_row(pivot).size() == 0,
        std::invalid_argument(
            "Chain_matrix::_remove_last - Column asked to be removed does not corresponds to a maximal simplex."));
    if constexpr (Master_matrix::Option_list::has_removable_rows) {
      RA_opt::erase_empty_row(pivot);
    }
  }
}
 
template <class Master_matrix>
inline void Chain_matrix<Master_matrix>::_update_barcode(Pos_index birth)
{
  if constexpr (Master_matrix::Option_list::has_column_pairings) {
    Pair_opt::_update_barcode(birth, nextPosition_);
  }
  ++nextPosition_;
}
 
template <class Master_matrix>
inline void Chain_matrix<Master_matrix>::_add_bar(Dimension dim)
{
  if constexpr (Master_matrix::Option_list::has_column_pairings) {
    Pair_opt::_add_bar(dim, nextPosition_);
  }
  ++nextPosition_;
}
 
template <class Master_matrix>
template <class Container>
inline void Chain_matrix<Master_matrix>::_container_insert(const Container& column, Index pos, Dimension dim)
{
  ID_index pivot = Master_matrix::get_row_index(*(column.rbegin()));
 
  if constexpr (Master_matrix::Option_list::has_map_column_container) {
    pivotToColumnIndex_.try_emplace(pivot, pos);
    if constexpr (Master_matrix::Option_list::has_row_access) {
      matrix_.try_emplace(pos, Column(pos, column, dim, RA_opt::_get_rows_ptr(), colSettings_));
    } else {
      matrix_.try_emplace(pos, Column(column, dim, colSettings_));
    }
  } else {
    if constexpr (Master_matrix::Option_list::has_row_access) {
      matrix_.emplace_back(pos, column, dim, RA_opt::_get_rows_ptr(), colSettings_);
    } else {
      matrix_.emplace_back(column, dim, colSettings_);
    }
    pivotToColumnIndex_[pivot] = pos;
  }
}
 
template <class Master_matrix>
template <class ColumnIterator> // Pair (pos,Column) if has_map_column_container, Column otherwise
inline void Chain_matrix<Master_matrix>::_container_insert(const ColumnIterator& rep)
{
  if constexpr (Master_matrix::Option_list::has_map_column_container) {
    const auto& col = rep.second;
    if constexpr (Master_matrix::Option_list::has_row_access) {
      matrix_.try_emplace(rep.first, Column(col, col.get_column_index(), RA_opt::_get_rows_ptr(), colSettings_));
    } else {
      matrix_.try_emplace(rep.first, Column(col, colSettings_));
    }
  } else {
    if constexpr (Master_matrix::Option_list::has_row_access) {
      matrix_.emplace_back(rep, rep.get_column_index(), RA_opt::_get_rows_ptr(), colSettings_);
    } else {
      matrix_.emplace_back(rep, colSettings_);
    }
  }
}
 
template <class Master_matrix>
inline void Chain_matrix<Master_matrix>::_insert_in(ID_index cellID, Tmp_column& column)
{
  if constexpr (Master_matrix::Option_list::is_z2)
    column.insert(cellID);
  else
    column.emplace(cellID, 1);
}
 
}  // namespace persistence_matrix
}  // namespace Gudhi
 
#endif  // PM_CHAIN_MATRIX_H