Boundary_matrix.h

Go to the documentation of this file.

/*    This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT.
 *    See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details.
 *    Author(s):       Hannah Schreiber
 *
 *    Copyright (C) 2022-24 Inria
 *
 *    Modification(s):
 *      - YYYY/MM Author: Description of the modification
 */
 
#ifndef PM_BOUNDARY_MATRIX_H
#define PM_BOUNDARY_MATRIX_H
 
#include <cassert>
#include <iostream> //print() only
#include <vector>
#include <utility>  //std::swap, std::move & std::exchange
 
namespace Gudhi {
namespace persistence_matrix {
 
// TODO: factorize/inherit/compose with Base_matrix?
template <class Master_matrix>
class Boundary_matrix : public Master_matrix::Matrix_dimension_option,
                        public Master_matrix::template Base_swap_option<Boundary_matrix<Master_matrix> >,
                        public Master_matrix::Base_pairing_option,
                        public Master_matrix::Matrix_row_access_option
{
 public:
  using Index = typename Master_matrix::Index;                        
  using ID_index = typename Master_matrix::ID_index;                  
  using Dimension = typename Master_matrix::Dimension;                
  using Field_operators = typename Master_matrix::Field_operators;
  using Field_element = typename Master_matrix::Element;                
  using Column = typename Master_matrix::Column;                        
  using Boundary = typename Master_matrix::Boundary;                    
  using Row = typename Master_matrix::Row;                              
  using Entry_constructor = typename Master_matrix::Entry_constructor;  
  using Column_settings = typename Master_matrix::Column_settings;      
  Boundary_matrix(Column_settings* colSettings);
  template <class Boundary_range = Boundary>
  Boundary_matrix(const std::vector<Boundary_range>& orderedBoundaries,
                  Column_settings* colSettings);
  Boundary_matrix(unsigned int numberOfColumns, Column_settings* colSettings);
  Boundary_matrix(const Boundary_matrix& matrixToCopy,
                  Column_settings* colSettings = nullptr);
  Boundary_matrix(Boundary_matrix&& other) noexcept;
 
  template <class Boundary_range = Boundary>
  Index insert_boundary(const Boundary_range& boundary,
                        Dimension dim = Master_matrix::template get_null_value<Dimension>());
  template <class Boundary_range = Boundary>
  Index insert_boundary(ID_index cellIndex, const Boundary_range& boundary,
                        Dimension dim = Master_matrix::template get_null_value<Dimension>());
  Column& get_column(Index columnIndex);
  Row& get_row(Index rowIndex);
  Index remove_last();
  void erase_empty_row(Index rowIndex);
 
  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);
  void zero_column(Index columnIndex);
  bool is_zero_entry(Index columnIndex, Index rowIndex) const;
  bool is_zero_column(Index columnIndex);
 
  Index get_pivot(Index columnIndex);
 
  void reset(Column_settings* colSettings) {
    matrix_.clear();
    nextInsertIndex_ = 0;
    colSettings_ = colSettings;
  }
 
  Boundary_matrix& operator=(const Boundary_matrix& other);
  friend void swap(Boundary_matrix& matrix1, Boundary_matrix& matrix2) {
    swap(static_cast<typename Master_matrix::Matrix_dimension_option&>(matrix1),
         static_cast<typename Master_matrix::Matrix_dimension_option&>(matrix2));
    swap(static_cast<typename Master_matrix::template Base_swap_option<Boundary_matrix<Master_matrix> >&>(matrix1),
         static_cast<typename Master_matrix::template Base_swap_option<Boundary_matrix<Master_matrix> >&>(matrix2));
    swap(static_cast<typename Master_matrix::Base_pairing_option&>(matrix1),
         static_cast<typename Master_matrix::Base_pairing_option&>(matrix2));
    matrix1.matrix_.swap(matrix2.matrix_);
    std::swap(matrix1.nextInsertIndex_, matrix2.nextInsertIndex_);
    std::swap(matrix1.colSettings_, matrix2.colSettings_);
 
    if constexpr (Master_matrix::Option_list::has_row_access) {
      swap(static_cast<typename Master_matrix::Matrix_row_access_option&>(matrix1),
           static_cast<typename Master_matrix::Matrix_row_access_option&>(matrix2));
    }
  }
 
  void print();  // for debug
 
 private:
  using Dim_opt = typename Master_matrix::Matrix_dimension_option;
  using Swap_opt = typename Master_matrix::template Base_swap_option<Boundary_matrix<Master_matrix> >;
  using Pair_opt = typename Master_matrix::Base_pairing_option;
  using RA_opt = typename Master_matrix::Matrix_row_access_option;
  using Column_container = typename Master_matrix::Column_container;
 
  friend Swap_opt;
  friend Pair_opt;
 
  Column_container matrix_;       
  Index nextInsertIndex_;         
  Column_settings* colSettings_;  
  static constexpr bool activeDimOption =
      Master_matrix::Option_list::has_matrix_maximal_dimension_access || Master_matrix::maxDimensionIsNeeded;
  static constexpr bool activeSwapOption =
      Master_matrix::Option_list::has_column_and_row_swaps || Master_matrix::Option_list::has_vine_update;
  static constexpr bool activePairingOption = Master_matrix::Option_list::has_column_pairings &&
                                              !Master_matrix::Option_list::has_vine_update &&
                                              !Master_matrix::Option_list::can_retrieve_representative_cycles;
 
  void _orderRowsIfNecessary();
  const Column& _get_column(Index columnIndex) const;
  Column& _get_column(Index columnIndex);
  Index _get_real_row_index(Index rowIndex) const;
  template <class Container>
  void _container_insert(const Container& column, Index pos, Dimension dim);
  void _container_insert(const Column& column, [[maybe_unused]] Index pos = 0);
};
 
template <class Master_matrix>
inline Boundary_matrix<Master_matrix>::Boundary_matrix(Column_settings* colSettings)
    : Dim_opt(Master_matrix::template get_null_value<Dimension>()),
      Swap_opt(),
      Pair_opt(),
      RA_opt(),
      nextInsertIndex_(0),
      colSettings_(colSettings)
{}
 
template <class Master_matrix>
template <class Boundary_range>
inline Boundary_matrix<Master_matrix>::Boundary_matrix(const std::vector<Boundary_range>& orderedBoundaries,
                                                       Column_settings* colSettings)
    : Dim_opt(Master_matrix::template get_null_value<Dimension>()),
      Swap_opt(orderedBoundaries.size()),
      Pair_opt(),
      RA_opt(orderedBoundaries.size()),
      nextInsertIndex_(orderedBoundaries.size()),
      colSettings_(colSettings)
{
  matrix_.reserve(orderedBoundaries.size());
 
  for (Index i = 0; i < orderedBoundaries.size(); i++) {
    _container_insert(orderedBoundaries[i], i, orderedBoundaries[i].size() == 0 ? 0 : orderedBoundaries[i].size() - 1);
  }
}
 
template <class Master_matrix>
inline Boundary_matrix<Master_matrix>::Boundary_matrix(unsigned int numberOfColumns,
                                                       Column_settings* colSettings)
    : Dim_opt(Master_matrix::template get_null_value<Dimension>()),
      Swap_opt(numberOfColumns),
      Pair_opt(),
      RA_opt(numberOfColumns),
      matrix_(!Master_matrix::Option_list::has_map_column_container && Master_matrix::Option_list::has_row_access
                  ? 0
                  : numberOfColumns),
      nextInsertIndex_(0),
      colSettings_(colSettings)
{
  if constexpr (!Master_matrix::Option_list::has_map_column_container && Master_matrix::Option_list::has_row_access)
    matrix_.reserve(numberOfColumns);
}
 
template <class Master_matrix>
inline Boundary_matrix<Master_matrix>::Boundary_matrix(const Boundary_matrix& matrixToCopy,
                                                       Column_settings* colSettings)
    : Dim_opt(static_cast<const Dim_opt&>(matrixToCopy)),
      Swap_opt(static_cast<const Swap_opt&>(matrixToCopy)),
      Pair_opt(static_cast<const Pair_opt&>(matrixToCopy)),
      RA_opt(static_cast<const RA_opt&>(matrixToCopy)),
      nextInsertIndex_(matrixToCopy.nextInsertIndex_),
      colSettings_(colSettings == nullptr ? matrixToCopy.colSettings_ : colSettings)
{
  matrix_.reserve(matrixToCopy.matrix_.size());
  for (const auto& cont : matrixToCopy.matrix_){
    if constexpr (Master_matrix::Option_list::has_map_column_container){
      _container_insert(cont.second, cont.first);
    } else {
      _container_insert(cont);
    }
  }
}
 
template <class Master_matrix>
inline Boundary_matrix<Master_matrix>::Boundary_matrix(Boundary_matrix&& other) noexcept
    : Dim_opt(std::move(static_cast<Dim_opt&>(other))),
      Swap_opt(std::move(static_cast<Swap_opt&>(other))),
      Pair_opt(std::move(static_cast<Pair_opt&>(other))),
      RA_opt(std::move(static_cast<RA_opt&>(other))),
      matrix_(std::move(other.matrix_)),
      nextInsertIndex_(std::exchange(other.nextInsertIndex_, 0)),
      colSettings_(std::exchange(other.colSettings_, nullptr))
{}
 
template <class Master_matrix>
template <class Boundary_range>
inline typename Boundary_matrix<Master_matrix>::Index Boundary_matrix<Master_matrix>::insert_boundary(
    const Boundary_range& boundary, Dimension dim)
{
  return insert_boundary(nextInsertIndex_, boundary, dim);
}
 
template <class Master_matrix>
template <class Boundary_range>
inline typename Boundary_matrix<Master_matrix>::Index Boundary_matrix<Master_matrix>::insert_boundary(
    ID_index cellIndex, const Boundary_range& boundary, Dimension dim)
{
  if (dim == Master_matrix::template get_null_value<Dimension>()) dim = boundary.size() == 0 ? 0 : boundary.size() - 1;
 
  _orderRowsIfNecessary();
 
  //updates container sizes
  if constexpr (Master_matrix::Option_list::has_row_access && !Master_matrix::Option_list::has_removable_rows) {
    if (boundary.size() != 0){
      ID_index pivot;
      if constexpr (Master_matrix::Option_list::is_z2) {
        pivot = *std::prev(boundary.end());
      } else {
        pivot = std::prev(boundary.end())->first;
      }
      //row container
      if (RA_opt::rows_->size() <= pivot) RA_opt::rows_->resize(pivot + 1);
    }
  }
 
  //row swap map containers
  if constexpr (Master_matrix::Option_list::has_map_column_container) {
    if constexpr (activeSwapOption) {
      Swap_opt::indexToRow_.emplace(cellIndex, cellIndex);
      Swap_opt::rowToIndex_.emplace(cellIndex, cellIndex);
    }
  } else {
    if constexpr (activeSwapOption) {
      for (Index i = Swap_opt::indexToRow_.size(); i <= cellIndex; ++i) {
        Swap_opt::indexToRow_.push_back(i);
        Swap_opt::rowToIndex_.push_back(i);
      }
    }
  }
 
  //maps for possible shifting between column content and position indices used for birth events
  if constexpr (activePairingOption){
    if (cellIndex != nextInsertIndex_){
      Pair_opt::idToPosition_.emplace(cellIndex, nextInsertIndex_);
      if constexpr (Master_matrix::Option_list::has_removable_columns){
        Pair_opt::PIDM::map_.emplace(nextInsertIndex_, cellIndex);
      }
    }
  }
 
  _container_insert(boundary, nextInsertIndex_, dim);
 
  return nextInsertIndex_++;
}
 
template <class Master_matrix>
inline typename Boundary_matrix<Master_matrix>::Column& Boundary_matrix<Master_matrix>::get_column(Index columnIndex)
{
  _orderRowsIfNecessary();
 
  return _get_column(columnIndex);
}
 
template <class Master_matrix>
inline typename Boundary_matrix<Master_matrix>::Row& Boundary_matrix<Master_matrix>::get_row(Index rowIndex)
{
  static_assert(Master_matrix::Option_list::has_row_access, "'get_row' is not implemented for the chosen options.");
 
  _orderRowsIfNecessary();
 
  return RA_opt::get_row(rowIndex);
}
 
template <class Master_matrix>
inline typename Boundary_matrix<Master_matrix>::Index Boundary_matrix<Master_matrix>::remove_last()
{
  static_assert(Master_matrix::Option_list::has_removable_columns,
                "'remove_last' is not implemented for the chosen options.");
 
  if (nextInsertIndex_ == 0) return Master_matrix::template get_null_value<Index>();  // empty matrix
  --nextInsertIndex_;
 
  //updates dimension max
  if constexpr (activeDimOption) {
    Dim_opt::update_down(matrix_.at(nextInsertIndex_).get_dimension());
  }
 
  //computes pivot and removes column from matrix_
  ID_index pivot;
  if constexpr (Master_matrix::Option_list::has_map_column_container) {
    auto it = matrix_.find(nextInsertIndex_);
    pivot = it->second.get_pivot();
    if constexpr (activeSwapOption) {
      // if the removed column is positive, the pivot won't change value
      if (Swap_opt::rowSwapped_ && pivot != Master_matrix::template get_null_value<ID_index>()) {
        Swap_opt::_orderRows();
        pivot = it->second.get_pivot();
      }
    }
    matrix_.erase(it);
  } else {
    pivot = matrix_[nextInsertIndex_].get_pivot();
    if constexpr (activeSwapOption) {
      // if the removed column is positive, the pivot won't change value
      if (Swap_opt::rowSwapped_ && pivot != Master_matrix::template get_null_value<ID_index>()) {
        Swap_opt::_orderRows();
        pivot = matrix_[nextInsertIndex_].get_pivot();
      }
    }
    if constexpr (Master_matrix::Option_list::has_row_access) {
      GUDHI_CHECK(nextInsertIndex_ == matrix_.size() - 1,
                  std::logic_error("Boundary_matrix::remove_last - Indexation problem."));
      matrix_.pop_back();
    } else {
      matrix_[nextInsertIndex_].clear();
    }
  }
 
  erase_empty_row(nextInsertIndex_);  // maximal, so empty
 
  //updates barcode
  if constexpr (activePairingOption) {
    Pair_opt::_remove_last(nextInsertIndex_);
  }
 
  return pivot;
}
 
template <class Master_matrix>
inline void Boundary_matrix<Master_matrix>::erase_empty_row(Index rowIndex)
{
  //computes real row index and erases it if necessary from the row swap map containers
  ID_index rowID = rowIndex;
  if constexpr (activeSwapOption) {
    if constexpr (Master_matrix::Option_list::has_map_column_container) {
      auto it = Swap_opt::indexToRow_.find(rowIndex);
      rowID = it->second;
      Swap_opt::rowToIndex_.erase(rowID);
      Swap_opt::indexToRow_.erase(it);
    } else {
      rowID = Swap_opt::indexToRow_[rowIndex];
    }
  }
 
  if constexpr (Master_matrix::Option_list::has_row_access && Master_matrix::Option_list::has_removable_rows) {
    RA_opt::erase_empty_row(rowID);
  }
}
 
template <class Master_matrix>
inline typename Boundary_matrix<Master_matrix>::Index Boundary_matrix<Master_matrix>::get_number_of_columns() const
{
  if constexpr (Master_matrix::Option_list::has_map_column_container) {
    return matrix_.size();
  } else {
    return nextInsertIndex_;  // matrix could have been resized much bigger while insert
  }
}
 
template <class Master_matrix>
inline typename Boundary_matrix<Master_matrix>::Dimension Boundary_matrix<Master_matrix>::get_column_dimension(
    Index columnIndex) const
{
  return _get_column(columnIndex).get_dimension();
}
 
template <class Master_matrix>
inline void Boundary_matrix<Master_matrix>::add_to(Index sourceColumnIndex, Index targetColumnIndex)
{
  _get_column(targetColumnIndex) += _get_column(sourceColumnIndex);
}
 
template <class Master_matrix>
inline void Boundary_matrix<Master_matrix>::multiply_target_and_add_to(Index sourceColumnIndex,
                                                                       const Field_element& coefficient,
                                                                       Index targetColumnIndex)
{
  _get_column(targetColumnIndex).multiply_target_and_add(coefficient, _get_column(sourceColumnIndex));
}
 
template <class Master_matrix>
inline void Boundary_matrix<Master_matrix>::multiply_source_and_add_to(const Field_element& coefficient,
                                                                       Index sourceColumnIndex,
                                                                       Index targetColumnIndex)
{
  _get_column(targetColumnIndex).multiply_source_and_add(_get_column(sourceColumnIndex), coefficient);
}
 
template <class Master_matrix>
inline void Boundary_matrix<Master_matrix>::zero_entry(Index columnIndex, Index rowIndex)
{
  _get_column(columnIndex).clear(_get_real_row_index(rowIndex));
}
 
template <class Master_matrix>
inline void Boundary_matrix<Master_matrix>::zero_column(Index columnIndex)
{
  _get_column(columnIndex).clear();
}
 
template <class Master_matrix>
inline bool Boundary_matrix<Master_matrix>::is_zero_entry(Index columnIndex, Index rowIndex) const
{
  return !(_get_column(columnIndex).is_non_zero(_get_real_row_index(rowIndex)));
}
 
template <class Master_matrix>
inline bool Boundary_matrix<Master_matrix>::is_zero_column(Index columnIndex)
{
  return _get_column(columnIndex).is_empty();
}
 
template <class Master_matrix>
inline typename Boundary_matrix<Master_matrix>::Index Boundary_matrix<Master_matrix>::get_pivot(Index columnIndex)
{
  _orderRowsIfNecessary();
 
  return _get_column(columnIndex).get_pivot();
}
 
template <class Master_matrix>
inline Boundary_matrix<Master_matrix>& Boundary_matrix<Master_matrix>::operator=(const Boundary_matrix& other)
{
  Dim_opt::operator=(other);
  Swap_opt::operator=(other);
  Pair_opt::operator=(other);
  RA_opt::operator=(other);
 
  matrix_.clear();
  nextInsertIndex_ = other.nextInsertIndex_;
  colSettings_ = other.colSettings_;
 
  matrix_.reserve(other.matrix_.size());
  for (const auto& cont : other.matrix_){
    if constexpr (Master_matrix::Option_list::has_map_column_container){
      _container_insert(cont.second, cont.first);
    } else {
      _container_insert(cont);
    }
  }
 
  return *this;
}
 
template <class Master_matrix>
inline void Boundary_matrix<Master_matrix>::print()
{
  if constexpr (activeSwapOption) {
    if (Swap_opt::rowSwapped_) Swap_opt::_orderRows();
  }
  std::cout << "Boundary_matrix:\n";
  for (Index i = 0; i < nextInsertIndex_; ++i) {
    Column& col = matrix_[i];
    for (auto e : col.get_content(nextInsertIndex_)) {
      if (e == 0u)
        std::cout << "- ";
      else
        std::cout << e << " ";
    }
    std::cout << "\n";
  }
  std::cout << "\n";
  if constexpr (Master_matrix::Option_list::has_row_access) {
    std::cout << "Row Matrix:\n";
    for (ID_index i = 0; i < nextInsertIndex_; ++i) {
      const auto& row = (*RA_opt::rows_)[i];
      for (const typename Column::Entry& entry : row) {
        std::cout << entry.get_column_index() << " ";
      }
      std::cout << "(" << i << ")\n";
    }
    std::cout << "\n";
  }
}
 
template <class Master_matrix>
inline void Boundary_matrix<Master_matrix>::_orderRowsIfNecessary()
{
  if constexpr (activeSwapOption) {
    if (Swap_opt::rowSwapped_) Swap_opt::_orderRows();
  }
}
 
template <class Master_matrix>
inline const typename Boundary_matrix<Master_matrix>::Column& Boundary_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 typename Boundary_matrix<Master_matrix>::Column& Boundary_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 typename Boundary_matrix<Master_matrix>::Index Boundary_matrix<Master_matrix>::_get_real_row_index(
    Index rowIndex) const
{
  if constexpr (Master_matrix::Option_list::has_column_and_row_swaps || Master_matrix::Option_list::has_vine_update) {
    if constexpr (Master_matrix::Option_list::has_map_column_container) {
      return Swap_opt::indexToRow_.at(rowIndex);
    } else {
      return Swap_opt::indexToRow_[rowIndex];
    }
  } else {
    return rowIndex;
  }
}
 
template <class Master_matrix>
template <class Container>
inline void Boundary_matrix<Master_matrix>::_container_insert(const Container& column,
                                                              Index pos,
                                                              Dimension dim)
{
  if constexpr (Master_matrix::Option_list::has_map_column_container) {
    if constexpr (Master_matrix::Option_list::has_row_access) {
      matrix_.try_emplace(pos, Column(pos, column, dim, RA_opt::rows_, 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::rows_, colSettings_);
    } else {
      if (matrix_.size() <= pos) {
        matrix_.emplace_back(column, dim, colSettings_);
      } else {
        matrix_[pos] = Column(column, dim, colSettings_);
      }
    }
  }
  if constexpr (activeDimOption) {
    Dim_opt::update_up(dim);
  }
}
 
template <class Master_matrix>
inline void Boundary_matrix<Master_matrix>::_container_insert(const Column& column, [[maybe_unused]] Index pos)
{
  if constexpr (Master_matrix::Option_list::has_map_column_container) {
    if constexpr (Master_matrix::Option_list::has_row_access) {
      matrix_.try_emplace(pos, Column(column, column.get_column_index(), RA_opt::rows_, colSettings_));
    } else {
      matrix_.try_emplace(pos, Column(column, colSettings_));
    }
  } else {
    if constexpr (Master_matrix::Option_list::has_row_access) {
      matrix_.emplace_back(column, column.get_column_index(), RA_opt::rows_, colSettings_);
    } else {
      matrix_.emplace_back(column, colSettings_);
    }
  }
}
 
}  // namespace persistence_matrix
}  // namespace Gudhi
 
#endif  // PM_BOUNDARY_MATRIX_H