RU_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 PM_RU_MATRIX_H
#define PM_RU_MATRIX_H
 
#include <utility>      //std::swap, std::move & std::exchange
#include <iostream>     //print() only
#include <vector>
#include <unordered_map>
 
namespace Gudhi {
namespace persistence_matrix {
 
template <class Master_matrix>
class RU_pairing;

template <class Master_matrix>
class RU_matrix : public Master_matrix::RU_pairing_option,
                  public Master_matrix::RU_vine_swap_option,
                  public Master_matrix::RU_representative_cycles_option
{
 private:
  using Pair_opt = typename Master_matrix::RU_pairing_option;
  using Swap_opt = typename Master_matrix::RU_vine_swap_option;
  using Rep_opt = typename Master_matrix::RU_representative_cycles_option;
 
 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_constructor = typename Master_matrix::Entry_constructor; 
  using Column_settings = typename Master_matrix::Column_settings;     
  using Boundary = typename Master_matrix::Boundary;                   
  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;                 

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

  template <class Boundary_range = Boundary>
  void insert_boundary(const Boundary_range& boundary,
                       Dimension dim = Master_matrix::template get_null_value<Dimension>());
  template <class Boundary_range = Boundary>
  void insert_boundary(ID_index cellIndex, const Boundary_range& boundary,
                       Dimension dim = Master_matrix::template get_null_value<Dimension>());
  template <class Boundary_range = Boundary>
  void insert_maximal_cell(Index columnIndex,
                           const Boundary_range& boundary,
                           Dimension dim = Master_matrix::template get_null_value<Dimension>());
  template <class Boundary_range = Boundary>
  void insert_maximal_cell(Index columnIndex,
                           ID_index cellIndex,
                           const Boundary_range& boundary,
                           Dimension dim = Master_matrix::template get_null_value<Dimension>());

  Column& get_column(Index columnIndex, bool inR = true);
  Row& get_row(Index rowIndex, bool inR = true);
  void erase_empty_row(Index rowIndex);
  void remove_maximal_cell(Index columnIndex);
  void remove_last();

  Dimension get_max_dimension() const;
  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);

  void zero_entry(Index columnIndex, Index rowIndex, bool inR = true);
  void zero_column(Index columnIndex, bool inR = true);
  bool is_zero_entry(Index columnIndex, Index rowIndex, bool inR = true) const;
  bool is_zero_column(Index columnIndex, bool inR = true);

  Index get_column_with_pivot(Index cellIndex) const;
  Index get_pivot(Index columnIndex);

  void reset(Column_settings* colSettings)
  {
    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();
    reducedMatrixR_.reset(colSettings);
    mirrorMatrixU_.reset(colSettings);
    pivotToColumnIndex_.clear();
    nextEventIndex_ = 0;
    positionToID_.clear();
    operators_ = Master_matrix::get_operator_ptr(colSettings);
  }

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

  friend void swap(RU_matrix& matrix1, RU_matrix& matrix2) noexcept
  {
    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(matrix1.reducedMatrixR_, matrix2.reducedMatrixR_);
    swap(matrix1.mirrorMatrixU_, matrix2.mirrorMatrixU_);
    matrix1.pivotToColumnIndex_.swap(matrix2.pivotToColumnIndex_);
    std::swap(matrix1.nextEventIndex_, matrix2.nextEventIndex_);
    matrix1.positionToID_.swap(matrix2.positionToID_);
    std::swap(matrix1.operators_, matrix2.operators_);
  }
 
  void print();  // for debug
 
 private:
  using Pivot_dictionary = typename Master_matrix::template Dictionary<Index>;
  using Position_dictionary = std::unordered_map<Pos_index, ID_index>;  // TODO: try other type of maps?
  using Barcode = typename Master_matrix::Barcode;
  using Bar_dictionary = typename Master_matrix::Bar_dictionary;
  using R_matrix = typename Master_matrix::Master_boundary_matrix;
  using U_matrix = typename Master_matrix::Master_base_matrix;
 
  friend Rep_opt;                    // direct access to the two matrices
  friend Swap_opt;                   // direct access to the two matrices, pivotToColumnIndex_
  friend RU_pairing<Master_matrix>;  // direct access to positionToID_
 
  R_matrix reducedMatrixR_; 
  // TODO: make U not accessible by default and add option to enable access? Inaccessible, it
  // needs less options and we could avoid some ifs.
  U_matrix mirrorMatrixU_;              
  Pivot_dictionary pivotToColumnIndex_; 
  Pos_index nextEventIndex_;            
  Position_dictionary positionToID_;    
  Field_operators const* operators_;    
 
  void _insert_boundary(Index currentIndex);
  void _initialize_U();
  void _reduce();
  void _reduce_last_column(Index lastIndex);
  void _reduce_column(Index target, Index eventIndex);
  void _reduce_column_by(Index target, Index source);
  Index _get_column_with_pivot(ID_index pivot) const;
  void _update_barcode(ID_index birthPivot, Pos_index death);
  void _add_bar(Dimension dim, Pos_index birth);
  void _remove_last_in_barcode(Pos_index eventIndex);
};
 
template <class Master_matrix>
inline RU_matrix<Master_matrix>::RU_matrix(Column_settings* colSettings)
    : Pair_opt(),
      Swap_opt(),
      Rep_opt(),
      reducedMatrixR_(colSettings),
      mirrorMatrixU_(colSettings),
      nextEventIndex_(0),
      operators_(Master_matrix::get_operator_ptr(colSettings))
{}
 
template <class Master_matrix>
template <class Boundary_range>
inline RU_matrix<Master_matrix>::RU_matrix(const std::vector<Boundary_range>& orderedBoundaries,
                                           Column_settings* colSettings)
    : Pair_opt(),
      Swap_opt(),
      Rep_opt(),
      reducedMatrixR_(orderedBoundaries, colSettings),
      mirrorMatrixU_(orderedBoundaries.size(), colSettings),
      nextEventIndex_(orderedBoundaries.size()),
      operators_(Master_matrix::get_operator_ptr(colSettings))
{
  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>());
  }
 
  _initialize_U();
  _reduce();
}
 
template <class Master_matrix>
inline RU_matrix<Master_matrix>::RU_matrix(unsigned int numberOfColumns, Column_settings* colSettings)
    : Pair_opt(),
      Swap_opt(),
      Rep_opt(),
      reducedMatrixR_(numberOfColumns, colSettings),
      mirrorMatrixU_(numberOfColumns, colSettings),
      nextEventIndex_(0),
      positionToID_(numberOfColumns),
      operators_(Master_matrix::get_operator_ptr(colSettings))
{
  if constexpr (Master_matrix::Option_list::has_map_column_container) {
    pivotToColumnIndex_.reserve(numberOfColumns);
  } else {
    pivotToColumnIndex_.resize(numberOfColumns, Master_matrix::template get_null_value<Index>());
  }
  if constexpr (Master_matrix::Option_list::has_column_pairings) {
    Pair_opt::_reserve(numberOfColumns);
  }
}
 
template <class Master_matrix>
inline RU_matrix<Master_matrix>::RU_matrix(const RU_matrix& matrixToCopy, Column_settings* colSettings)
    : Pair_opt(static_cast<const Pair_opt&>(matrixToCopy)),
      Swap_opt(static_cast<const Swap_opt&>(matrixToCopy)),
      Rep_opt(static_cast<const Rep_opt&>(matrixToCopy)),
      reducedMatrixR_(matrixToCopy.reducedMatrixR_, colSettings),
      mirrorMatrixU_(matrixToCopy.mirrorMatrixU_, colSettings),
      pivotToColumnIndex_(matrixToCopy.pivotToColumnIndex_),
      nextEventIndex_(matrixToCopy.nextEventIndex_),
      positionToID_(matrixToCopy.positionToID_),
      operators_(colSettings == nullptr ? matrixToCopy.operators_ : Master_matrix::get_operator_ptr(colSettings))
{}
 
template <class Master_matrix>
inline RU_matrix<Master_matrix>::RU_matrix(RU_matrix&& other) noexcept
    : 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))),
      reducedMatrixR_(std::move(other.reducedMatrixR_)),
      mirrorMatrixU_(std::move(other.mirrorMatrixU_)),
      pivotToColumnIndex_(std::move(other.pivotToColumnIndex_)),
      nextEventIndex_(std::exchange(other.nextEventIndex_, 0)),
      positionToID_(std::move(other.positionToID_)),
      operators_(std::exchange(other.operators_, nullptr))
{}
 
template <class Master_matrix>
template <class Boundary_range>
inline void RU_matrix<Master_matrix>::insert_boundary(const Boundary_range& boundary, Dimension dim)
{
  _insert_boundary(reducedMatrixR_.insert_boundary(boundary, dim));
}
 
template <class Master_matrix>
template <class Boundary_range>
inline void RU_matrix<Master_matrix>::insert_boundary(ID_index cellIndex, const Boundary_range& boundary, Dimension dim)
{
  // maps for possible shifting between column content and position indices used for birth events
  if (cellIndex != nextEventIndex_) {
    positionToID_.emplace(nextEventIndex_, cellIndex);
    if constexpr (Master_matrix::Option_list::has_column_pairings) {
      Pair_opt::_insert_id_position(cellIndex, nextEventIndex_);
    }
  }
 
  _insert_boundary(reducedMatrixR_.insert_boundary(cellIndex, boundary, dim));
}
 
template <class Master_matrix>
template <class Boundary_range>
inline void RU_matrix<Master_matrix>::insert_maximal_cell(Index columnIndex,
                                                          const Boundary_range& boundary,
                                                          Dimension dim)
{
  static_assert(Master_matrix::Option_list::has_vine_update,
                "'insert_maximal_cell' is not implemented for the chosen options.");
 
  GUDHI_CHECK(columnIndex >= 0, std::invalid_argument("Indices have to be positive."));
 
  insert_boundary(boundary, dim);
 
  // If started with 0 columns, no swaps are needed
  if (get_number_of_columns() == 1) return;
 
  for (Index curr = get_number_of_columns() - 1; curr > columnIndex; --curr) {
    Swap_opt::vine_swap(curr - 1);
  }
}
 
template <class Master_matrix>
template <class Boundary_range>
inline void RU_matrix<Master_matrix>::insert_maximal_cell(Index columnIndex,
                                                          ID_index cellIndex,
                                                          const Boundary_range& boundary,
                                                          Dimension dim)
{
  static_assert(Master_matrix::Option_list::has_vine_update,
                "'insert_maximal_cell' is not implemented for the chosen options.");
 
  GUDHI_CHECK(columnIndex >= 0, std::invalid_argument("Indices have to be positive."));
 
  insert_boundary(cellIndex, boundary, dim);
 
  // TODO: factorize following with the other insert_maximal_cell ?
 
  // If started with 0 columns, no swaps are needed
  if (get_number_of_columns() == 1) return;
 
  for (Index curr = get_number_of_columns() - 1; curr > columnIndex; --curr) {
    Swap_opt::vine_swap(curr - 1);
  }
}
 
template <class Master_matrix>
inline typename RU_matrix<Master_matrix>::Column& RU_matrix<Master_matrix>::get_column(Index columnIndex, bool inR)
{
  if (inR) {
    return reducedMatrixR_.get_column(columnIndex);
  }
  return mirrorMatrixU_.get_column(columnIndex);
}
 
template <class Master_matrix>
inline typename RU_matrix<Master_matrix>::Row& RU_matrix<Master_matrix>::get_row(Index rowIndex, bool inR)
{
  static_assert(Master_matrix::Option_list::has_row_access, "'get_row' is not implemented for the chosen options.");
 
  if (inR) {
    return reducedMatrixR_.get_row(rowIndex);
  }
  return mirrorMatrixU_.get_row(rowIndex);
}
 
template <class Master_matrix>
inline void RU_matrix<Master_matrix>::erase_empty_row(Index rowIndex)
{
  reducedMatrixR_.erase_empty_row(rowIndex);
}
 
template <class Master_matrix>
inline void RU_matrix<Master_matrix>::remove_maximal_cell(Index columnIndex)
{
  static_assert(Master_matrix::Option_list::has_removable_columns && Master_matrix::Option_list::has_vine_update,
                "'remove_maximal_cell' is not implemented for the chosen options.");
 
  // TODO: is there an easy test to verify maximality even without row access?
 
  for (Index curr = columnIndex; curr < nextEventIndex_ - 1; ++curr) {
    Swap_opt::vine_swap(curr);
  }
 
  remove_last();
}
 
template <class Master_matrix>
inline void RU_matrix<Master_matrix>::remove_last()
{
  static_assert(Master_matrix::Option_list::has_removable_columns,
                "'remove_last' is not implemented for the chosen options.");
 
  if (nextEventIndex_ == 0) return;  // empty matrix
  --nextEventIndex_;
 
  // assumes PosIdx == MatIdx for boundary matrices.
  _remove_last_in_barcode(nextEventIndex_);
 
  mirrorMatrixU_.remove_last();
  if constexpr (Master_matrix::Option_list::has_map_column_container) {
    pivotToColumnIndex_.erase(reducedMatrixR_.remove_last());
  } else {
    ID_index lastPivot = reducedMatrixR_.remove_last();
    if (lastPivot != Master_matrix::template get_null_value<ID_index>())
      pivotToColumnIndex_[lastPivot] = Master_matrix::template get_null_value<Index>();
  }
 
  // if has_column_pairings is true, then the element is already removed in _remove_last_in_barcode
  // to avoid a second "find"
  if constexpr (!Master_matrix::Option_list::has_column_pairings) {
    positionToID_.erase(nextEventIndex_);
  }
}
 
template <class Master_matrix>
inline typename RU_matrix<Master_matrix>::Dimension RU_matrix<Master_matrix>::get_max_dimension() const
{
  return reducedMatrixR_.get_max_dimension();
}
 
template <class Master_matrix>
inline typename RU_matrix<Master_matrix>::Index RU_matrix<Master_matrix>::get_number_of_columns() const
{
  return reducedMatrixR_.get_number_of_columns();
}
 
template <class Master_matrix>
inline typename RU_matrix<Master_matrix>::Dimension RU_matrix<Master_matrix>::get_column_dimension(
    Index columnIndex) const
{
  return reducedMatrixR_.get_column_dimension(columnIndex);
}
 
template <class Master_matrix>
inline void RU_matrix<Master_matrix>::add_to(Index sourceColumnIndex, Index targetColumnIndex)
{
  reducedMatrixR_.add_to(sourceColumnIndex, targetColumnIndex);
  // U transposed to avoid row operations
  if constexpr (Master_matrix::Option_list::is_z2)
    mirrorMatrixU_.add_to(targetColumnIndex, sourceColumnIndex);
  else
    mirrorMatrixU_.multiply_source_and_add_to(
        operators_->get_characteristic() - 1, targetColumnIndex, sourceColumnIndex);
}
 
template <class Master_matrix>
inline void RU_matrix<Master_matrix>::multiply_target_and_add_to(Index sourceColumnIndex,
                                                                 const Field_element& coefficient,
                                                                 Index targetColumnIndex)
{
  reducedMatrixR_.multiply_target_and_add_to(sourceColumnIndex, coefficient, targetColumnIndex);
  // U transposed to avoid row operations
  mirrorMatrixU_.get_column(targetColumnIndex) *= coefficient;
  mirrorMatrixU_.multiply_source_and_add_to(operators_->get_characteristic() - 1, targetColumnIndex, sourceColumnIndex);
}
 
template <class Master_matrix>
inline void RU_matrix<Master_matrix>::multiply_source_and_add_to(const Field_element& coefficient,
                                                                 Index sourceColumnIndex,
                                                                 Index targetColumnIndex)
{
  reducedMatrixR_.multiply_source_and_add_to(coefficient, sourceColumnIndex, targetColumnIndex);
  // U transposed to avoid row operations
  if constexpr (Master_matrix::Option_list::is_z2) {
    if (coefficient) mirrorMatrixU_.add_to(targetColumnIndex, sourceColumnIndex);
  } else {
    mirrorMatrixU_.multiply_source_and_add_to(
        operators_->get_characteristic() - coefficient, targetColumnIndex, sourceColumnIndex);
  }
}
 
template <class Master_matrix>
inline void RU_matrix<Master_matrix>::zero_entry(Index columnIndex, Index rowIndex, bool inR)
{
  if (inR) {
    return reducedMatrixR_.zero_entry(columnIndex, rowIndex);
  }
  return mirrorMatrixU_.zero_entry(columnIndex, rowIndex);
}
 
template <class Master_matrix>
inline void RU_matrix<Master_matrix>::zero_column(Index columnIndex, bool inR)
{
  if (inR) {
    return reducedMatrixR_.zero_column(columnIndex);
  }
  return mirrorMatrixU_.zero_column(columnIndex);
}
 
template <class Master_matrix>
inline bool RU_matrix<Master_matrix>::is_zero_entry(Index columnIndex, Index rowIndex, bool inR) const
{
  if (inR) {
    return reducedMatrixR_.is_zero_entry(columnIndex, rowIndex);
  }
  return mirrorMatrixU_.is_zero_entry(columnIndex, rowIndex);
}
 
template <class Master_matrix>
inline bool RU_matrix<Master_matrix>::is_zero_column(Index columnIndex, bool inR)
{
  if (inR) {
    return reducedMatrixR_.is_zero_column(columnIndex);
  }
  return mirrorMatrixU_.is_zero_column(columnIndex);
}
 
template <class Master_matrix>
inline typename RU_matrix<Master_matrix>::Index RU_matrix<Master_matrix>::get_column_with_pivot(Index cellIndex) const
{
  if constexpr (Master_matrix::Option_list::has_map_column_container) {
    return pivotToColumnIndex_.at(cellIndex);
  } else {
    return pivotToColumnIndex_[cellIndex];
  }
}
 
template <class Master_matrix>
inline typename RU_matrix<Master_matrix>::Index RU_matrix<Master_matrix>::get_pivot(Index columnIndex)
{
  return reducedMatrixR_.get_column(columnIndex).get_pivot();
}
 
template <class Master_matrix>
inline RU_matrix<Master_matrix>& RU_matrix<Master_matrix>::operator=(const RU_matrix& other)
{
  if (this == &other) return *this;
 
  Swap_opt::operator=(other);
  Pair_opt::operator=(other);
  Rep_opt::operator=(other);
  reducedMatrixR_ = other.reducedMatrixR_;
  mirrorMatrixU_ = other.mirrorMatrixU_;
  pivotToColumnIndex_ = other.pivotToColumnIndex_;
  nextEventIndex_ = other.nextEventIndex_;
  positionToID_ = other.positionToID_;
  operators_ = other.operators_;
 
  return *this;
}
 
template <class Master_matrix>
inline RU_matrix<Master_matrix>& RU_matrix<Master_matrix>::operator=(RU_matrix&& other) noexcept
{
  if (this == &other) return *this;
 
  Pair_opt::operator=(std::move(other));
  Swap_opt::operator=(std::move(other));
  Rep_opt::operator=(std::move(other));
 
  reducedMatrixR_ = std::move(other.reducedMatrixR_);
  mirrorMatrixU_ = std::move(other.mirrorMatrixU_);
  pivotToColumnIndex_ = std::move(other.pivotToColumnIndex_);
  nextEventIndex_ = std::exchange(other.nextEventIndex_, 0);
  positionToID_ = std::move(other.positionToID_);
  operators_ = std::exchange(other.operators_, nullptr);
}
 
template <class Master_matrix>
inline void RU_matrix<Master_matrix>::print()
{
  std::cout << "R_matrix:\n";
  reducedMatrixR_.print();
  std::cout << "U_matrix:\n";
  mirrorMatrixU_.print();
}
 
template <class Master_matrix>
inline void RU_matrix<Master_matrix>::_insert_boundary(Index currentIndex)
{
  mirrorMatrixU_.insert_column(currentIndex, 1);
 
  if constexpr (!Master_matrix::Option_list::has_map_column_container) {
    ID_index pivot = reducedMatrixR_.get_column(currentIndex).get_pivot();
    if (pivot != Master_matrix::template get_null_value<ID_index>() && pivotToColumnIndex_.size() <= pivot)
      pivotToColumnIndex_.resize((pivot + 1) * 2, Master_matrix::template get_null_value<Index>());
  }
 
  _reduce_last_column(currentIndex);
  ++nextEventIndex_;
}
 
template <class Master_matrix>
inline void RU_matrix<Master_matrix>::_initialize_U()
{
  for (ID_index i = 0; i < reducedMatrixR_.get_number_of_columns(); i++) {
    mirrorMatrixU_.insert_column(i, 1);
  }
}
 
template <class Master_matrix>
inline void RU_matrix<Master_matrix>::_reduce()
{
  if constexpr (Master_matrix::Option_list::has_column_pairings) {
    Pair_opt::_reserve(reducedMatrixR_.get_number_of_columns());
  }
 
  for (Index i = 0; i < reducedMatrixR_.get_number_of_columns(); i++) {
    if (!(reducedMatrixR_.is_zero_column(i))) {
      _reduce_column(i, i);
    } else {
      _add_bar(get_column_dimension(i), i);
    }
  }
}
 
template <class Master_matrix>
inline void RU_matrix<Master_matrix>::_reduce_last_column(Index lastIndex)
{
  if (reducedMatrixR_.get_column(lastIndex).is_empty()) {
    _add_bar(get_column_dimension(lastIndex), nextEventIndex_);
    return;
  }
 
  _reduce_column(lastIndex, nextEventIndex_);
}
 
template <class Master_matrix>
inline void RU_matrix<Master_matrix>::_reduce_column(Index target, Index eventIndex)
{
  Column& curr = reducedMatrixR_.get_column(target);
  ID_index pivot = curr.get_pivot();
  Index currIndex = _get_column_with_pivot(pivot);
 
  while (pivot != Master_matrix::template get_null_value<ID_index>() &&
         currIndex != Master_matrix::template get_null_value<Index>()) {
    _reduce_column_by(target, currIndex);
    pivot = curr.get_pivot();
    currIndex = _get_column_with_pivot(pivot);
  }
 
  if (pivot != Master_matrix::template get_null_value<ID_index>()) {
    if constexpr (Master_matrix::Option_list::has_map_column_container) {
      pivotToColumnIndex_.try_emplace(pivot, target);
    } else {
      pivotToColumnIndex_[pivot] = target;
    }
    _update_barcode(pivot, eventIndex);
  } else {
    _add_bar(get_column_dimension(target), eventIndex);
  }
}
 
template <class Master_matrix>
inline void RU_matrix<Master_matrix>::_reduce_column_by(Index target, Index source)
{
  Column& curr = reducedMatrixR_.get_column(target);
  if constexpr (Master_matrix::Option_list::is_z2) {
    curr += reducedMatrixR_.get_column(source);
    // to avoid having to do line operations, U is transposed
    // TODO: explain this somewhere in the documentation...
    mirrorMatrixU_.get_column(source).push_back(*mirrorMatrixU_.get_column(target).begin());
  } else {
    Column& toadd = reducedMatrixR_.get_column(source);
    Field_element coef = toadd.get_pivot_value();
    coef = operators_->get_inverse(coef);
    operators_->multiply_inplace(coef, operators_->get_characteristic() - curr.get_pivot_value());
 
    curr.multiply_source_and_add(toadd, coef);
    auto entry = *mirrorMatrixU_.get_column(target).begin();
    operators_->multiply_inplace(entry.get_element(), operators_->get_characteristic() - coef);
    // to avoid having to do line operations, U is transposed
    // TODO: explain this somewhere in the documentation...
    mirrorMatrixU_.get_column(source).push_back(entry);
  }
}
 
template <class Master_matrix>
inline typename RU_matrix<Master_matrix>::Index RU_matrix<Master_matrix>::_get_column_with_pivot(ID_index pivot) const
{
  if (pivot == Master_matrix::template get_null_value<ID_index>())
    return Master_matrix::template get_null_value<Index>();
  if constexpr (Master_matrix::Option_list::has_map_column_container) {
    auto it = pivotToColumnIndex_.find(pivot);
    if (it == pivotToColumnIndex_.end()) return Master_matrix::template get_null_value<Index>();
    return it->second;
  } else {
    if (pivot >= pivotToColumnIndex_.size()) return Master_matrix::template get_null_value<Index>();
    return pivotToColumnIndex_[pivot];
  }
}
 
template <class Master_matrix>
inline void RU_matrix<Master_matrix>::_update_barcode(ID_index birthPivot, Pos_index death)
{
  if constexpr (Master_matrix::Option_list::has_column_pairings) {
    Pair_opt::_update_barcode(birthPivot, death);
  }
}
 
template <class Master_matrix>
inline void RU_matrix<Master_matrix>::_add_bar(Dimension dim, Pos_index birth)
{
  if constexpr (Master_matrix::Option_list::has_column_pairings) {
    Pair_opt::_add_bar(dim, birth);
  }
}
 
template <class Master_matrix>
inline void RU_matrix<Master_matrix>::_remove_last_in_barcode(Pos_index eventIndex)
{
  if constexpr (Master_matrix::Option_list::has_column_pairings) {
    Pair_opt::_remove_last(eventIndex);
  }
}
 
}  // namespace persistence_matrix
}  // namespace Gudhi
 
#endif  // PM_RU_MATRIX_H