Base_matrix_with_column_compression.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_BASE_MATRIX_COMPRESSION_H
#define PM_BASE_MATRIX_COMPRESSION_H
 
#include <iostream>  //print() only
#include <vector>
#include <utility>  //std::swap, std::move & std::exchange
 
#include <boost/intrusive/set.hpp>
#include <boost/pending/disjoint_sets.hpp>
 
#include <gudhi/Simple_object_pool.h>
 
namespace Gudhi {
namespace persistence_matrix {

template <class Master_matrix>
class Base_matrix_with_column_compression : protected Master_matrix::Matrix_row_access_option
{
 public:
  using Index = typename Master_matrix::Index;         
  using Dimension = typename Master_matrix::Dimension; 
  using Field_operators = typename Master_matrix::Field_operators;
  using Field_element = typename Master_matrix::Element;               
  using Row = typename Master_matrix::Row;                             
  using Entry_constructor = typename Master_matrix::Entry_constructor; 
  using Column_settings = typename Master_matrix::Column_settings;     

  class Column : public Master_matrix::Column,
                 public boost::intrusive::set_base_hook<boost::intrusive::link_mode<boost::intrusive::normal_link> >
  {
   public:
    using Base = typename Master_matrix::Column;
 
    Column(Column_settings* colSettings = nullptr) : Base(colSettings) {}
 
    template <class Container>
    Column(const Container& nonZeroRowIndices, Column_settings* colSettings) : Base(nonZeroRowIndices, colSettings)
    {}
 
    template <class Container, class Row_container>
    Column(Index columnIndex,
           const Container& nonZeroRowIndices,
           Row_container* rowContainer,
           Column_settings* colSettings)
        : Base(columnIndex, nonZeroRowIndices, rowContainer, colSettings)
    {}
 
    template <class Container, class = std::enable_if_t<!std::is_arithmetic_v<Container> > >
    Column(const Container& nonZeroRowIndices, Dimension dimension, Column_settings* colSettings)
        : Base(nonZeroRowIndices, dimension, colSettings)
    {}
 
    template <class Container, class Row_container, class = std::enable_if_t<!std::is_arithmetic_v<Container> > >
    Column(Index columnIndex,
           const Container& nonZeroRowIndices,
           Dimension dimension,
           Row_container* rowContainer,
           Column_settings* colSettings)
        : Base(columnIndex, nonZeroRowIndices, dimension, rowContainer, colSettings)
    {}
 
    Column(Index idx, Dimension dimension, Column_settings* colSettings) : Base(idx, dimension, colSettings) {}
 
    Column(Index idx, Field_element e, Dimension dimension, Column_settings* colSettings)
        : Base(idx, e, dimension, colSettings)
    {}
 
    template <class Row_container>
    Column(Index columnIndex, Index idx, Dimension dimension, Row_container* rowContainer, Column_settings* colSettings)
        : Base(columnIndex, idx, dimension, rowContainer, colSettings)
    {}
 
    template <class Row_container>
    Column(Index columnIndex,
           Index idx,
           Field_element e,
           Dimension dimension,
           Row_container* rowContainer,
           Column_settings* colSettings)
        : Base(columnIndex, idx, e, dimension, rowContainer, colSettings)
    {}
 
    Column(const Column& column, Column_settings* colSettings = nullptr)
        : Base(static_cast<const Base&>(column), colSettings), rep_(column.rep_)
    {}
 
    template <class Row_container>
    Column(const Column& column, Index columnIndex, Row_container* rowContainer, Column_settings* colSettings = nullptr)
        : Base(static_cast<const Base&>(column), columnIndex, rowContainer, colSettings), rep_(column.rep_)
    {}
 
    Column(Column&& column) noexcept : Base(std::move(static_cast<Base&>(column))), rep_(std::exchange(column.rep_, 0))
    {}
 
    ~Column() = default;
 
    // TODO: is it possible to make this work?
    //  template <class... U>
    //  Column(U&&... u) : Base(std::forward<U>(u)...) {}
 
    Index get_rep() const { return rep_; }
 
    void set_rep(const Index& rep) { rep_ = rep; }
 
    Column& operator=(const Column& other)
    {
      Base::operator=(other);
      rep_ = other.rep_;
    }
 
    Column& operator=(Column&& other) noexcept
    {
      Base::operator=(std::move(other));
      rep_ = std::exchange(other.rep_, 0);
    }
 
    struct Hasher {
      size_t operator()(const Column& c) const { return std::hash<Base>()(static_cast<Base>(c)); }
    };
 
   private:
    Index rep_; 
  };

  Base_matrix_with_column_compression(Column_settings* colSettings);
  template <class Container>
  Base_matrix_with_column_compression(const std::vector<Container>& columns, Column_settings* colSettings);
  Base_matrix_with_column_compression(unsigned int numberOfColumns, Column_settings* colSettings);
  Base_matrix_with_column_compression(const Base_matrix_with_column_compression& matrixToCopy,
                                      Column_settings* colSettings = nullptr);
  Base_matrix_with_column_compression(Base_matrix_with_column_compression&& other) noexcept;
  ~Base_matrix_with_column_compression();

  template <class Container, class = std::enable_if_t<!std::is_arithmetic_v<Container> > >
  void insert_column(const Container& column);
  void insert_column(Index idx);
  void insert_column(Index idx, Field_element e);
  template <class Boundary_range>
  void insert_boundary(const Boundary_range& boundary,
                       Dimension dim = Master_matrix::template get_null_value<Dimension>());
  const Column& get_column(Index columnIndex);
  const Row& get_row(Index rowIndex) const;
  void erase_empty_row(Index rowIndex);

  Index get_number_of_columns() const;

  template <class Entry_range_or_column_index>
  void add_to(const Entry_range_or_column_index& sourceColumn, Index targetColumnIndex);
  template <class Entry_range_or_column_index>
  void multiply_target_and_add_to(const Entry_range_or_column_index& sourceColumn,
                                  const Field_element& coefficient,
                                  Index targetColumnIndex);
  template <class Entry_range_or_column_index>
  void multiply_source_and_add_to(const Field_element& coefficient,
                                  const Entry_range_or_column_index& sourceColumn,
                                  Index targetColumnIndex);

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

  void reset(Column_settings* colSettings)
  {
    columnToRep_.clear_and_dispose(Delete_disposer(this));
    columnClasses_ = boost::disjoint_sets_with_storage<>();
    repToColumn_.clear();
    nextColumnIndex_ = 0;
    colSettings_ = colSettings;
  }

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

  friend void swap(Base_matrix_with_column_compression& matrix1, Base_matrix_with_column_compression& matrix2) noexcept
  {
    matrix1.columnToRep_.swap(matrix2.columnToRep_);
    std::swap(matrix1.columnClasses_, matrix2.columnClasses_);
    matrix1.repToColumn_.swap(matrix2.repToColumn_);
    std::swap(matrix1.nextColumnIndex_, matrix2.nextColumnIndex_);
    std::swap(matrix1.colSettings_, matrix2.colSettings_);
    std::swap(matrix1.columnPool_, matrix2.columnPool_);
 
    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:
  struct Delete_disposer {
    Delete_disposer(Base_matrix_with_column_compression* matrix) : matrix_(matrix) {}
 
    void operator()(Column* delete_this) { matrix_->columnPool_->destroy(delete_this); }
 
    Base_matrix_with_column_compression* matrix_;
  };
 
  using RA_opt = typename Master_matrix::Matrix_row_access_option;
  using Col_dict = boost::intrusive::set<Column, boost::intrusive::constant_time_size<false> >;
 
  Col_dict columnToRep_;                              
  boost::disjoint_sets_with_storage<> columnClasses_; 
  std::vector<Column*> repToColumn_;                  
  Index nextColumnIndex_;                             
  Column_settings* colSettings_;                      
  std::unique_ptr<Simple_object_pool<Column> > columnPool_;
  inline static const Column empty_column_; 
 
  template <class F>
  void _initialize_column(F&& get_column);
  void _insert_column(Index columnIndex);
  void _insert_double_column(Index columnIndex, typename Col_dict::iterator& doubleIt);
};
 
template <class Master_matrix>
inline Base_matrix_with_column_compression<Master_matrix>::Base_matrix_with_column_compression(
    Column_settings* colSettings)
    : RA_opt(), nextColumnIndex_(0), colSettings_(colSettings), columnPool_(new Simple_object_pool<Column>)
{
}
 
template <class Master_matrix>
template <class Container>
inline Base_matrix_with_column_compression<Master_matrix>::Base_matrix_with_column_compression(
    const std::vector<Container>& columns,
    Column_settings* colSettings)
    : RA_opt(columns.size()),
      columnClasses_(columns.size()),
      repToColumn_(columns.size(), nullptr),
      nextColumnIndex_(0),
      colSettings_(colSettings),
      columnPool_(new Simple_object_pool<Column>)
{
  for (const Container& c : columns) {
    insert_column(c);
  }
}
 
template <class Master_matrix>
inline Base_matrix_with_column_compression<Master_matrix>::Base_matrix_with_column_compression(
    unsigned int numberOfColumns,
    Column_settings* colSettings)
    : RA_opt(numberOfColumns),
      columnClasses_(numberOfColumns),
      repToColumn_(numberOfColumns, nullptr),
      nextColumnIndex_(0),
      colSettings_(colSettings),
      columnPool_(new Simple_object_pool<Column>)
{
}
 
template <class Master_matrix>
inline Base_matrix_with_column_compression<Master_matrix>::Base_matrix_with_column_compression(
    const Base_matrix_with_column_compression& matrixToCopy,
    Column_settings* colSettings)
    : RA_opt(static_cast<const RA_opt&>(matrixToCopy)),
      columnClasses_(matrixToCopy.columnClasses_),
      repToColumn_(matrixToCopy.repToColumn_.size(), nullptr),
      nextColumnIndex_(0),
      colSettings_(colSettings == nullptr ? matrixToCopy.colSettings_ : colSettings),
      columnPool_(new Simple_object_pool<Column>)
{
  for (const Column* col : matrixToCopy.repToColumn_) {
    if (col != nullptr) {
      if constexpr (Master_matrix::Option_list::has_row_access) {
        repToColumn_[nextColumnIndex_] =
            columnPool_->construct(*col, col->get_column_index(), RA_opt::_get_rows_ptr(), colSettings_);
      } else {
        repToColumn_[nextColumnIndex_] = columnPool_->construct(*col, colSettings_);
      }
      columnToRep_.insert(columnToRep_.end(), *repToColumn_[nextColumnIndex_]);
      repToColumn_[nextColumnIndex_]->set_rep(nextColumnIndex_);
    }
    nextColumnIndex_++;
  }
}
 
template <class Master_matrix>
inline Base_matrix_with_column_compression<Master_matrix>::Base_matrix_with_column_compression(
    Base_matrix_with_column_compression&& other) noexcept
    : RA_opt(std::move(static_cast<RA_opt&>(other))),
      columnToRep_(std::move(other.columnToRep_)),
      columnClasses_(std::move(other.columnClasses_)),
      repToColumn_(std::move(other.repToColumn_)),
      nextColumnIndex_(std::exchange(other.nextColumnIndex_, 0)),
      colSettings_(std::exchange(other.colSettings_, nullptr)),
      columnPool_(std::exchange(other.columnPool_, nullptr))
{
}
 
template <class Master_matrix>
inline Base_matrix_with_column_compression<Master_matrix>::~Base_matrix_with_column_compression()
{
  columnToRep_.clear_and_dispose(Delete_disposer(this));
}
 
template <class Master_matrix>
template <class Container, class>
inline void Base_matrix_with_column_compression<Master_matrix>::insert_column(const Container& column)
{
  insert_boundary(column);
}
 
template <class Master_matrix>
inline void Base_matrix_with_column_compression<Master_matrix>::insert_column(Index idx)
{
  static_assert(Master_matrix::Option_list::is_z2,
                "Insertion method not available for Zp != Z2. Please specify the coefficient.");
 
  if constexpr (Master_matrix::Option_list::has_row_access && !Master_matrix::Option_list::has_removable_rows) {
    RA_opt::_resize(idx);
  }
 
  _initialize_column([&]() -> Column* {
    if constexpr (Master_matrix::Option_list::has_row_access) {
      return columnPool_->construct(nextColumnIndex_, idx, 0, RA_opt::_get_rows_ptr(), colSettings_);
    } else {
      return columnPool_->construct(idx, 0, colSettings_);
    }
  });
 
  _insert_column(nextColumnIndex_);
 
  nextColumnIndex_++;
}
 
template <class Master_matrix>
inline void Base_matrix_with_column_compression<Master_matrix>::insert_column(Index idx,
                                                                              Field_element e)
{
  static_assert(!Master_matrix::Option_list::is_z2,
                "Insertion method not available for Zp == Z2. Please do not specify any coefficient.");
 
  if constexpr (Master_matrix::Option_list::has_row_access && !Master_matrix::Option_list::has_removable_rows) {
    RA_opt::_resize(idx);
  }
 
  _initialize_column([&]() -> Column* {
    if constexpr (Master_matrix::Option_list::has_row_access) {
      return columnPool_->construct(nextColumnIndex_, idx, e, 0, RA_opt::_get_rows_ptr(), colSettings_);
    } else {
      return columnPool_->construct(idx, e, 0, colSettings_);
    }
  });
 
  _insert_column(nextColumnIndex_);
 
  nextColumnIndex_++;
}
 
template <class Master_matrix>
template <class Boundary_range>
inline void Base_matrix_with_column_compression<Master_matrix>::insert_boundary(const Boundary_range& boundary,
                                                                                Dimension dim)
{
  // handles a dimension which is not actually stored.
  // TODO: verify if this is not a problem for the cohomology algorithm and if yes,
  // change how Column_dimension_option is chosen.
  if (dim == Master_matrix::template get_null_value<Dimension>()) dim = boundary.size() == 0 ? 0 : boundary.size() - 1;
 
  if constexpr (Master_matrix::Option_list::has_row_access && !Master_matrix::Option_list::has_removable_rows) {
    if (boundary.begin() != boundary.end()) {
      RA_opt::_resize(Master_matrix::get_row_index(*std::prev(boundary.end())));
    }
  }
 
  _initialize_column([&]() -> Column* {
    if constexpr (Master_matrix::Option_list::has_row_access) {
      return columnPool_->construct(nextColumnIndex_, boundary, dim, RA_opt::_get_rows_ptr(), colSettings_);
    } else {
      return columnPool_->construct(boundary, dim, colSettings_);
    }
  });
 
  _insert_column(nextColumnIndex_);
 
  nextColumnIndex_++;
}
 
template <class Master_matrix>
inline const typename Base_matrix_with_column_compression<Master_matrix>::Column&
Base_matrix_with_column_compression<Master_matrix>::get_column(Index columnIndex)
{
  auto col = repToColumn_[columnClasses_.find_set(columnIndex)];
  if (col == nullptr) return empty_column_;
  return *col;
}
 
template <class Master_matrix>
inline const typename Base_matrix_with_column_compression<Master_matrix>::Row&
Base_matrix_with_column_compression<Master_matrix>::get_row(Index rowIndex) const
{
  static_assert(Master_matrix::Option_list::has_row_access, "Row access has to be enabled for this method.");
 
  return RA_opt::get_row(rowIndex);
}
 
template <class Master_matrix>
inline void Base_matrix_with_column_compression<Master_matrix>::erase_empty_row(Index rowIndex)
{
  if constexpr (Master_matrix::Option_list::has_row_access && Master_matrix::Option_list::has_removable_rows) {
    RA_opt::erase_empty_row(rowIndex);
  }
}
 
template <class Master_matrix>
inline typename Base_matrix_with_column_compression<Master_matrix>::Index
Base_matrix_with_column_compression<Master_matrix>::get_number_of_columns() const
{
  return nextColumnIndex_;
}
 
template <class Master_matrix>
template <class Entry_range_or_column_index>
inline void Base_matrix_with_column_compression<Master_matrix>::add_to(const Entry_range_or_column_index& sourceColumn,
                                                                       Index targetColumnIndex)
{
  // handle case where targetRep == sourceRep?
  Index targetRep = columnClasses_.find_set(targetColumnIndex);
  Column& target = *repToColumn_[targetRep];
  columnToRep_.erase(target);
  if constexpr (std::is_integral_v<Entry_range_or_column_index>) {
    target += get_column(sourceColumn);
  } else {
    target += sourceColumn;
  }
  _insert_column(targetRep);
}
 
template <class Master_matrix>
template <class Entry_range_or_column_index>
inline void Base_matrix_with_column_compression<Master_matrix>::multiply_target_and_add_to(
    const Entry_range_or_column_index& sourceColumn,
    const Field_element& coefficient,
    Index targetColumnIndex)
{
  // handle case where targetRep == sourceRep?
  Index targetRep = columnClasses_.find_set(targetColumnIndex);
  Column& target = *repToColumn_[targetRep];
  columnToRep_.erase(target);
  if constexpr (std::is_integral_v<Entry_range_or_column_index>) {
    target.multiply_target_and_add(coefficient, get_column(sourceColumn));
  } else {
    target.multiply_target_and_add(coefficient, sourceColumn);
  }
  _insert_column(targetRep);
}
 
template <class Master_matrix>
template <class Entry_range_or_column_index>
inline void Base_matrix_with_column_compression<Master_matrix>::multiply_source_and_add_to(
    const Field_element& coefficient,
    const Entry_range_or_column_index& sourceColumn,
    Index targetColumnIndex)
{
  // handle case where targetRep == sourceRep?
  Index targetRep = columnClasses_.find_set(targetColumnIndex);
  Column& target = *repToColumn_[targetRep];
  columnToRep_.erase(target);
  if constexpr (std::is_integral_v<Entry_range_or_column_index>) {
    target.multiply_source_and_add(get_column(sourceColumn), coefficient);
  } else {
    target.multiply_source_and_add(sourceColumn, coefficient);
  }
  _insert_column(targetRep);
}
 
template <class Master_matrix>
inline bool Base_matrix_with_column_compression<Master_matrix>::is_zero_entry(Index columnIndex, Index rowIndex)
{
  auto col = repToColumn_[columnClasses_.find_set(columnIndex)];
  if (col == nullptr) return true;
  return !col->is_non_zero(rowIndex);
}
 
template <class Master_matrix>
inline bool Base_matrix_with_column_compression<Master_matrix>::is_zero_column(Index columnIndex)
{
  auto col = repToColumn_[columnClasses_.find_set(columnIndex)];
  if (col == nullptr) return true;
  return col->is_empty();
}
 
template <class Master_matrix>
inline Base_matrix_with_column_compression<Master_matrix>&
Base_matrix_with_column_compression<Master_matrix>::operator=(const Base_matrix_with_column_compression& other)
{
  if (this == &other) return *this;
 
  for (auto col : repToColumn_) {
    if (col != nullptr) {
      columnPool_->destroy(col);
      col = nullptr;
    }
  }
  RA_opt::operator=(other);
  columnClasses_ = other.columnClasses_;
  columnToRep_.reserve(other.columnToRep_.size());
  repToColumn_.resize(other.repToColumn_.size(), nullptr);
  nextColumnIndex_ = 0;
  colSettings_ = other.colSettings_;
  for (const Column* col : other.repToColumn_) {
    if constexpr (Master_matrix::Option_list::has_row_access) {
      repToColumn_[nextColumnIndex_] =
          columnPool_->construct(*col, col->get_column_index(), RA_opt::_get_rows_ptr(), colSettings_);
    } else {
      repToColumn_[nextColumnIndex_] = columnPool_->construct(*col, colSettings_);
    }
    columnToRep_.insert(columnToRep_.end(), *repToColumn_[nextColumnIndex_]);
    repToColumn_[nextColumnIndex_]->set_rep(nextColumnIndex_);
    nextColumnIndex_++;
  }
  return *this;
}
 
template <class Master_matrix>
inline Base_matrix_with_column_compression<Master_matrix>&
Base_matrix_with_column_compression<Master_matrix>::operator=(Base_matrix_with_column_compression&& other) noexcept
{
  if (&repToColumn_ == &(other.repToColumn_)) return *this;
 
  RA_opt::operator=(std::move(other));
 
  columnToRep_.clear_and_dispose(Delete_disposer(this));
 
  columnToRep_ = std::move(other.columnToRep_);
  columnClasses_ = std::move(other.columnClasses_);
  repToColumn_ = std::move(other.repToColumn_);
  nextColumnIndex_ = std::exchange(other.nextColumnIndex_, 0);
  colSettings_ = std::exchange(other.colSettings_, nullptr);
  columnPool_ = std::exchange(other.columnPool_, nullptr);
 
  return *this;
}
 
template <class Master_matrix>
inline void Base_matrix_with_column_compression<Master_matrix>::print()
{
  std::cout << "Compressed_matrix:\n";
  for (Column& col : columnToRep_) {
    for (auto e : col->get_content(nextColumnIndex_)) {
      if (e == 0U)
        std::cout << "- ";
      else
        std::cout << e << " ";
    }
    std::cout << "(";
    for (Index i = 0; i < nextColumnIndex_; ++i) {
      if (columnClasses_.find_set(i) == col.get_rep()) std::cout << i << " ";
    }
    std::cout << ")\n";
  }
  std::cout << "\n";
  std::cout << "Row Matrix:\n";
  for (Index i = 0; i < RA_opt::_get_rows_ptr()->size(); ++i) {
    const Row& row = RA_opt::get_row(i);
    for (const auto& entry : row) {
      std::cout << entry.get_column_index() << " ";
    }
    std::cout << "(" << i << ")\n";
  }
  std::cout << "\n";
}
 
template <class Master_matrix>
template <class F>
inline void Base_matrix_with_column_compression<Master_matrix>::_initialize_column(F&& get_column)
{
  if (repToColumn_.size() == nextColumnIndex_) {
    // could perhaps be avoided, if find_set returns something special when it does not find
    columnClasses_.link(nextColumnIndex_, nextColumnIndex_);
    repToColumn_.push_back(std::forward<F>(get_column)());
  } else {
    repToColumn_[nextColumnIndex_] = std::forward<F>(get_column)();
  }
}
 
template <class Master_matrix>
inline void Base_matrix_with_column_compression<Master_matrix>::_insert_column(Index columnIndex)
{
  Column& col = *repToColumn_[columnIndex];
 
  if (col.is_empty()) {
    columnPool_->destroy(&col);  // delete curr_col;
    repToColumn_[columnIndex] = nullptr;
    return;
  }
 
  col.set_rep(columnIndex);
  auto res = columnToRep_.insert(col);
  if (res.first->get_rep() != columnIndex) {  // if true, then redundant column
    _insert_double_column(columnIndex, res.first);
  }
}
 
template <class Master_matrix>
inline void Base_matrix_with_column_compression<Master_matrix>::_insert_double_column(
    Index columnIndex,
    typename Col_dict::iterator& doubleIt)
{
  Index doubleRep = doubleIt->get_rep();
  columnClasses_.link(columnIndex, doubleRep);  // both should be representatives
  Index newRep = columnClasses_.find_set(columnIndex);
 
  columnPool_->destroy(repToColumn_[columnIndex]);
  repToColumn_[columnIndex] = nullptr;
 
  if (newRep == columnIndex) {
    std::swap(repToColumn_[doubleRep], repToColumn_[columnIndex]);
    doubleIt->set_rep(columnIndex);
  }
}
 
}  // namespace persistence_matrix
}  // namespace Gudhi
 
#endif  // PM_BASE_MATRIX_COMPRESSION_H