Merge branch 'master' into pip_package_for_python_3_13

CMakeLists.txt

@@ -47,6 +47,7 @@ add_gudhi_module(Toplex_map)
 add_gudhi_module(Witness_complex)
 add_gudhi_module(Nerve_GIC)
 add_gudhi_module(Persistence_matrix)
+add_gudhi_module(Zigzag_persistence)
 
 # Include module CMake subdirectories
 # GUDHI_SUB_DIRECTORIES is managed in CMAKE_MODULE_PATH/GUDHI_modules.cmake

azure-pipelines.yml

@@ -4,7 +4,7 @@ jobs:
     timeoutInMinutes: 0
     cancelTimeoutInMinutes: 60
     pool:
-      vmImage: macOS-latest
+      vmImage: macOS-13
     variables:
       pythonVersion: "3.9"
       cmakeBuildType: Release

src/CMakeLists.txt

@@ -46,6 +46,7 @@ add_gudhi_module(Toplex_map)
 add_gudhi_module(Witness_complex)
 add_gudhi_module(Nerve_GIC)
 add_gudhi_module(Persistence_matrix)
+add_gudhi_module(Zigzag_persistence)
 
 set(GUDHI_BIBLIO_DIR ${CMAKE_SOURCE_DIR})
 # For "make doxygen" - Requires GUDHI_USER_VERSION_DIR to be set

src/Persistence_matrix/concept/PersistenceMatrixColumn.h

@@ -426,6 +426,15 @@ class PersistenceMatrixColumn :
   template <class Entry_range>
   PersistenceMatrixColumn& multiply_source_and_add(const Entry_range& column, const Field_element& val);
 
+  /**
+   * @brief Adds a copy of the given entry at the end of the column. It is therefore assumed that the row index
+   * of the entry is higher than the current pivot of the column. Not available for @ref chainmatrix "chain matrices"
+   * and is only needed for @ref rumatrix "RU matrices".
+   * 
+   * @param entry Entry to push back.
+   */
+  void push_back(const Entry& entry);
+
   /**
    * @brief Equality comparator. Equal in the sense that what is "supposed" to be contained in the columns is equal,
    * not what is actually stored in the underlying container. For example, the underlying container of 

src/Persistence_matrix/include/gudhi/Persistence_matrix/Base_matrix.h

@@ -615,7 +615,7 @@ inline void Base_matrix<Master_matrix>::print()
   if constexpr (Master_matrix::Option_list::has_row_access) {
     std::cout << "Row Matrix:\n";
     for (Index i = 0; i < nextInsertIndex_; ++i) {
-      const auto& row = RA_opt::rows_[i];
+      const auto& row = (*RA_opt::rows_)[i];
       for (const auto& entry : row) {
         std::cout << entry.get_column_index() << " ";
       }

src/Persistence_matrix/include/gudhi/Persistence_matrix/Boundary_matrix.h

@@ -729,8 +729,8 @@ inline void Boundary_matrix<Master_matrix>::print()
   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 auto& entry : row) {
+      const auto& row = (*RA_opt::rows_)[i];
+      for (const typename Column::Entry& entry : row) {
         std::cout << entry.get_column_index() << " ";
       }
       std::cout << "(" << i << ")\n";

src/Persistence_matrix/include/gudhi/Persistence_matrix/Chain_matrix.h

@@ -913,24 +913,28 @@ 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() && pivotToColumnIndex_[i] != static_cast<Index>(-1); ++i) {
+    for (ID_index i = 0; i < pivotToColumnIndex_.size(); ++i) {
       Index pos = pivotToColumnIndex_[i];
-      const Column& col = matrix_[pos];
-      for (const auto& entry : col) {
-        std::cout << entry.get_row_index() << " ";
+      if (pos != static_cast<Index>(-1)){
+        const Column& col = matrix_[pos];
+        for (const auto& entry : col) {
+          std::cout << entry.get_row_index() << " ";
+        }
+        std::cout << "(" << i << ", " << pos << ")\n";
       }
-      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() && pivotToColumnIndex_[i] != static_cast<Index>(-1); ++i) {
+      for (ID_index i = 0; i < pivotToColumnIndex_.size(); ++i) {
         Index pos = pivotToColumnIndex_[i];
-        const Row& row = RA_opt::get_row(pos);
-        for (const auto& entry : row) {
-          std::cout << entry.get_column_index() << " ";
+        if (pos != static_cast<Index>(-1)){
+          const Row& row = RA_opt::get_row(pos);
+          for (const auto& entry : row) {
+            std::cout << entry.get_column_index() << " ";
+          }
+          std::cout << "(" << i << ", " << pos << ")\n";
         }
-        std::cout << "(" << i << ", " << pos << ")\n";
       }
     }
   } else {

src/Persistence_matrix/include/gudhi/Persistence_matrix/RU_matrix.h

@@ -538,7 +538,7 @@ inline void RU_matrix<Master_matrix>::insert_boundary(ID_index cellIndex,
   if constexpr (Master_matrix::Option_list::has_vine_update) {
     if (cellIndex != nextEventIndex_) {
       Swap_opt::_positionToRowIdx().emplace(nextEventIndex_, cellIndex);
-      if (Master_matrix::Option_list::has_column_pairings) {
+      if constexpr (Master_matrix::Option_list::has_column_pairings) {
         Swap_opt::template RU_pairing<Master_matrix>::idToPosition_.emplace(cellIndex, nextEventIndex_);
       }
     }
@@ -638,7 +638,7 @@ inline void RU_matrix<Master_matrix>::add_to(Index sourceColumnIndex, Index targ
 {
   reducedMatrixR_.add_to(sourceColumnIndex, targetColumnIndex);
   // U transposed to avoid row operations
-  if constexpr (Master_matrix::Option_list::has_vine_update)
+  if constexpr (Master_matrix::Option_list::is_z2)
     mirrorMatrixU_.add_to(targetColumnIndex, sourceColumnIndex);
   else
     mirrorMatrixU_.add_to(sourceColumnIndex, targetColumnIndex);
@@ -649,6 +649,8 @@ inline void RU_matrix<Master_matrix>::multiply_target_and_add_to(Index sourceCol
                                                                  const Field_element& coefficient,
                                                                  Index targetColumnIndex)
 {
+  static_assert(!Master_matrix::Option_list::is_z2,
+                "Multiplication with something else than the identity is not allowed with Z2 coefficients.");
   reducedMatrixR_.multiply_target_and_add_to(sourceColumnIndex, coefficient, targetColumnIndex);
   mirrorMatrixU_.multiply_target_and_add_to(sourceColumnIndex, coefficient, targetColumnIndex);
 }
@@ -658,6 +660,8 @@ inline void RU_matrix<Master_matrix>::multiply_source_and_add_to(const Field_ele
                                                                  Index sourceColumnIndex,
                                                                  Index targetColumnIndex)
 {
+  static_assert(!Master_matrix::Option_list::is_z2,
+                "Multiplication with something else than the identity is not allowed with Z2 coefficients.");
   reducedMatrixR_.multiply_source_and_add_to(coefficient, sourceColumnIndex, targetColumnIndex);
   mirrorMatrixU_.multiply_source_and_add_to(coefficient, sourceColumnIndex, targetColumnIndex);
 }
@@ -844,19 +848,17 @@ inline void RU_matrix<Master_matrix>::_reduce_column_by(Index target, Index sour
     curr += reducedMatrixR_.get_column(source);
     // to avoid having to do line operations during vineyards, U is transposed
     // TODO: explain this somewhere in the documentation...
-    if constexpr (Master_matrix::Option_list::has_vine_update)
-      mirrorMatrixU_.get_column(source) += mirrorMatrixU_.get_column(target);
-    else
-      mirrorMatrixU_.get_column(target) += mirrorMatrixU_.get_column(source);
+    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);
+    // but no transposition for Zp, careful if there will be vineyard or rep cycles in Zp one day
+    // TODO: explain this somewhere in the documentation...
     mirrorMatrixU_.multiply_source_and_add_to(coef, source, target);
-    // mirrorMatrixU_.get_column(target).multiply_source_and_add(mirrorMatrixU_.get_column(source), coef);
   }
 }
 

src/Persistence_matrix/include/gudhi/Persistence_matrix/base_pairing.h

@@ -132,10 +132,11 @@ inline const typename Base_pairing<Master_matrix>::Barcode& Base_pairing<Master_
 template <class Master_matrix>
 inline void Base_pairing<Master_matrix>::_reduce() 
 {
-  std::unordered_map<ID_index, Index> pivotsToColumn(_matrix()->get_number_of_columns());
+  std::unordered_map<Index, Index> negativeColumns(_matrix()->get_number_of_columns());
 
   auto dim = _matrix()->get_max_dimension();
   std::vector<std::vector<Index> > columnsByDim(dim + 1);
+  for (auto& v : columnsByDim) v.reserve(_matrix()->get_number_of_columns());
   for (unsigned int i = 0; i < _matrix()->get_number_of_columns(); i++) {
     columnsByDim[dim - _matrix()->get_column_dimension(i)].push_back(i);
   }
@@ -144,17 +145,21 @@ inline void Base_pairing<Master_matrix>::_reduce()
     for (Index i : cols) {
       auto& curr = _matrix()->get_column(i);
       if (curr.is_empty()) {
-        if (pivotsToColumn.find(i) == pivotsToColumn.end()) {
+        if (negativeColumns.find(i) == negativeColumns.end()) {
           barcode_.emplace_back(i, -1, dim);
         }
       } else {
         ID_index pivot = curr.get_pivot();
+        auto it = idToPosition_.find(pivot);
+        Index pivotColumnNumber = it == idToPosition_.end() ? pivot : it->second;
+        auto itNeg = negativeColumns.find(pivotColumnNumber);
+        Index pivotKiller = itNeg == negativeColumns.end() ? -1 : itNeg->second;
 
-        while (pivot != static_cast<ID_index>(-1) && pivotsToColumn.find(pivot) != pivotsToColumn.end()) {
+        while (pivot != static_cast<ID_index>(-1) && pivotKiller != static_cast<Index>(-1)) {
           if constexpr (Master_matrix::Option_list::is_z2) {
-            curr += _matrix()->get_column(pivotsToColumn.at(pivot));
+            curr += _matrix()->get_column(pivotKiller);
           } else {
-            auto& toadd = _matrix()->get_column(pivotsToColumn.at(pivot));
+            auto& toadd = _matrix()->get_column(pivotKiller);
             typename Master_matrix::Element coef = toadd.get_pivot_value();
             auto& operators = _matrix()->colSettings_->operators;
             coef = operators.get_inverse(coef);
@@ -163,12 +168,14 @@ inline void Base_pairing<Master_matrix>::_reduce()
           }
 
           pivot = curr.get_pivot();
+          it = idToPosition_.find(pivot);
+          pivotColumnNumber = it == idToPosition_.end() ? pivot : it->second;
+          itNeg = negativeColumns.find(pivotColumnNumber);
+          pivotKiller = itNeg == negativeColumns.end() ? -1 : itNeg->second;
         }
 
         if (pivot != static_cast<ID_index>(-1)) {
-          pivotsToColumn.emplace(pivot, i);
-          auto it = idToPosition_.find(pivot);
-          auto pivotColumnNumber = it == idToPosition_.end() ? pivot : it->second;
+          negativeColumns.emplace(pivotColumnNumber, i);
           _matrix()->get_column(pivotColumnNumber).clear();
           barcode_.emplace_back(pivotColumnNumber, i, dim - 1);
         } else {

src/Persistence_matrix/include/gudhi/Persistence_matrix/columns/heap_column.h

@@ -25,6 +25,7 @@
 #include <utility>    //std::swap, std::move & std::exchange
 
 #include <boost/iterator/indirect_iterator.hpp>
+#include "gudhi/Debug_utils.h"
 
 #include <gudhi/Persistence_matrix/allocators/entry_constructors.h>
 
@@ -134,6 +135,8 @@ class Heap_column : public Master_matrix::Column_dimension_option, public Master
   Heap_column& multiply_source_and_add(const Entry_range& column, const Field_element& val);
   Heap_column& multiply_source_and_add(Heap_column& column, const Field_element& val);
 
+  void push_back(const Entry& entry);
+
   std::size_t compute_hash_value();
 
   friend bool operator==(const Heap_column& c1, const Heap_column& c2) {
@@ -868,6 +871,21 @@ inline Heap_column<Master_matrix>& Heap_column<Master_matrix>::multiply_source_a
   return *this;
 }
 
+template <class Master_matrix>
+inline void Heap_column<Master_matrix>::push_back(const Entry& entry)
+{
+  static_assert(Master_matrix::Option_list::is_of_boundary_type, "`push_back` is not available for Chain matrices.");
+
+  GUDHI_CHECK(entry.get_row_index() > get_pivot(), "The new row index has to be higher than the current pivot.");
+
+  Entry* newEntry = entryPool_->construct(entry.get_row_index());
+  if constexpr (!Master_matrix::Option_list::is_z2) {
+    newEntry->set_element(operators_->get_value(entry.get_element()));
+  }
+  column_.push_back(newEntry);
+  std::push_heap(column_.begin(), column_.end(), entryPointerComp_);
+}
+
 template <class Master_matrix>
 inline Heap_column<Master_matrix>& Heap_column<Master_matrix>::operator=(const Heap_column& other)
 {

src/Persistence_matrix/include/gudhi/Persistence_matrix/columns/intrusive_list_column.h

@@ -132,6 +132,8 @@ class Intrusive_list_column : public Master_matrix::Row_access_option,
   Intrusive_list_column& multiply_source_and_add(const Entry_range& column, const Field_element& val);
   Intrusive_list_column& multiply_source_and_add(Intrusive_list_column& column, const Field_element& val);
 
+  void push_back(const Entry& entry);
+
   friend bool operator==(const Intrusive_list_column& c1, const Intrusive_list_column& c2) {
     if (&c1 == &c2) return true;
 
@@ -821,6 +823,20 @@ inline Intrusive_list_column<Master_matrix>& Intrusive_list_column<Master_matrix
   return *this;
 }
 
+template <class Master_matrix>
+inline void Intrusive_list_column<Master_matrix>::push_back(const Entry& entry)
+{
+  static_assert(Master_matrix::Option_list::is_of_boundary_type, "`push_back` is not available for Chain matrices.");
+
+  GUDHI_CHECK(entry.get_row_index() > get_pivot(), "The new row index has to be higher than the current pivot.");
+
+  if constexpr (Master_matrix::Option_list::is_z2) {
+    _insert_entry(entry.get_row_index(), column_.end());
+  } else {
+    _insert_entry(entry.get_element(), entry.get_row_index(), column_.end());
+  }
+}
+
 template <class Master_matrix>
 inline Intrusive_list_column<Master_matrix>& Intrusive_list_column<Master_matrix>::operator=(
     const Intrusive_list_column& other)

src/Persistence_matrix/include/gudhi/Persistence_matrix/columns/intrusive_set_column.h

@@ -133,6 +133,8 @@ class Intrusive_set_column : public Master_matrix::Row_access_option,
   Intrusive_set_column& multiply_source_and_add(const Entry_range& column, const Field_element& val);
   Intrusive_set_column& multiply_source_and_add(Intrusive_set_column& column, const Field_element& val);
 
+  void push_back(const Entry& entry);
+
   friend bool operator==(const Intrusive_set_column& c1, const Intrusive_set_column& c2) {
     if (&c1 == &c2) return true;
 
@@ -821,6 +823,20 @@ inline Intrusive_set_column<Master_matrix>& Intrusive_set_column<Master_matrix>:
   return *this;
 }
 
+template <class Master_matrix>
+inline void Intrusive_set_column<Master_matrix>::push_back(const Entry& entry)
+{
+  static_assert(Master_matrix::Option_list::is_of_boundary_type, "`push_back` is not available for Chain matrices.");
+
+  GUDHI_CHECK(entry.get_row_index() > get_pivot(), "The new row index has to be higher than the current pivot.");
+
+  if constexpr (Master_matrix::Option_list::is_z2) {
+    _insert_entry(entry.get_row_index(), column_.end());
+  } else {
+    _insert_entry(entry.get_element(), entry.get_row_index(), column_.end());
+  }
+}
+
 template <class Master_matrix>
 inline Intrusive_set_column<Master_matrix>& Intrusive_set_column<Master_matrix>::operator=(
     const Intrusive_set_column& other)

src/Persistence_matrix/include/gudhi/Persistence_matrix/columns/list_column.h

@@ -131,6 +131,8 @@ class List_column : public Master_matrix::Row_access_option,
   List_column& multiply_source_and_add(const Entry_range& column, const Field_element& val);
   List_column& multiply_source_and_add(List_column& column, const Field_element& val);
 
+  void push_back(const Entry& entry);
+
   friend bool operator==(const List_column& c1, const List_column& c2) {
     if (&c1 == &c2) return true;
 
@@ -805,6 +807,20 @@ inline List_column<Master_matrix>& List_column<Master_matrix>::multiply_source_a
   return *this;
 }
 
+template <class Master_matrix>
+inline void List_column<Master_matrix>::push_back(const Entry& entry)
+{
+  static_assert(Master_matrix::Option_list::is_of_boundary_type, "`push_back` is not available for Chain matrices.");
+
+  GUDHI_CHECK(entry.get_row_index() > get_pivot(), "The new row index has to be higher than the current pivot.");
+
+  if constexpr (Master_matrix::Option_list::is_z2) {
+    _insert_entry(entry.get_row_index(), column_.end());
+  } else {
+    _insert_entry(entry.get_element(), entry.get_row_index(), column_.end());
+  }
+}
+
 template <class Master_matrix>
 inline List_column<Master_matrix>& List_column<Master_matrix>::operator=(const List_column& other)
 {

src/Persistence_matrix/include/gudhi/Persistence_matrix/columns/naive_vector_column.h

@@ -131,6 +131,8 @@ class Naive_vector_column : public Master_matrix::Row_access_option,
   Naive_vector_column& multiply_source_and_add(const Entry_range& column, const Field_element& val);
   Naive_vector_column& multiply_source_and_add(Naive_vector_column& column, const Field_element& val);
 
+  void push_back(const Entry& entry);
+
   friend bool operator==(const Naive_vector_column& c1, const Naive_vector_column& c2) {
     if (&c1 == &c2) return true;
     if (c1.column_.size() != c2.column_.size()) return false;
@@ -801,6 +803,20 @@ inline Naive_vector_column<Master_matrix>& Naive_vector_column<Master_matrix>::m
   return *this;
 }
 
+template <class Master_matrix>
+inline void Naive_vector_column<Master_matrix>::push_back(const Entry& entry)
+{
+  static_assert(Master_matrix::Option_list::is_of_boundary_type, "`push_back` is not available for Chain matrices.");
+
+  GUDHI_CHECK(entry.get_row_index() > get_pivot(), "The new row index has to be higher than the current pivot.");
+
+  if constexpr (Master_matrix::Option_list::is_z2) {
+    _insert_entry(entry.get_row_index(), column_);
+  } else {
+    _insert_entry(entry.get_element(), entry.get_row_index(), column_);
+  }
+}
+
 template <class Master_matrix>
 inline Naive_vector_column<Master_matrix>& Naive_vector_column<Master_matrix>::operator=(
     const Naive_vector_column& other)