Merge pull request GUDHI#1109 from hschreiber/matrix_barcode_fix

[Persistence matrix] barcode index shift fix

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

@@ -23,6 +23,8 @@
 #include <algorithm>
 #include <vector>
 
+#include "boundary_face_position_to_id_mapper.h"
+
 namespace Gudhi {
 namespace persistence_matrix {
 
@@ -46,7 +48,11 @@ struct Dummy_base_pairing {
  * @tparam Master_matrix An instantiation of @ref Matrix from which all types and options are deduced.
  */
 template <class Master_matrix>
-class Base_pairing 
+class Base_pairing : public std::conditional<
+                       Master_matrix::Option_list::has_removable_columns,
+                       Face_position_to_ID_mapper<typename Master_matrix::id_index, typename Master_matrix::pos_index>,
+                       Dummy_pos_mapper
+                    >::type
 {
  public:
   using Bar = typename Master_matrix::Bar;                            /**< Bar type. */
@@ -59,18 +65,6 @@ class Base_pairing
    * @brief Default constructor.
    */
   Base_pairing();
-  /**
-   * @brief Copy constructor.
-   * 
-   * @param matrixToCopy Matrix to copy.
-   */
-  Base_pairing(const Base_pairing& matrixToCopy);
-  /**
-   * @brief Move constructor.
-   * 
-   * @param other Matrix to move.
-   */
-  Base_pairing(Base_pairing&& other) noexcept;
 
   /**
    * @brief Reduces the matrix stored in @ref Boundary_matrix and computes the corresponding barcode.
@@ -83,26 +77,37 @@ class Base_pairing
    */
   const barcode_type& get_current_barcode();
 
-  /**
-   * @brief Assign operator.
-   */
-  Base_pairing& operator=(Base_pairing other);
   /**
    * @brief Swap operator.
    */
   friend void swap(Base_pairing& pairing1, Base_pairing& pairing2) {
+    if constexpr (Master_matrix::Option_list::has_removable_columns) {
+      swap(static_cast<Face_position_to_ID_mapper<id_index, pos_index>&>(pairing1),
+           static_cast<Face_position_to_ID_mapper<id_index, pos_index>&>(pairing2));
+    }
     pairing1.barcode_.swap(pairing2.barcode_);
     pairing1.deathToBar_.swap(pairing2.deathToBar_);
+    pairing1.idToPosition_.swap(pairing2.idToPosition_);
     std::swap(pairing1.isReduced_, pairing2.isReduced_);
   }
 
  protected:
   using pos_index = typename Master_matrix::pos_index;
+  using id_index = typename Master_matrix::id_index;
   using dictionary_type = typename Master_matrix::bar_dictionary_type;
   using base_matrix = typename Master_matrix::Boundary_matrix_type;
+  //PIDM = Position to ID Map
+  using PIDM = typename std::conditional<Master_matrix::Option_list::has_removable_columns,
+                                        Face_position_to_ID_mapper<id_index, pos_index>,
+                                        Dummy_pos_mapper
+                                       >::type;
 
   barcode_type barcode_;        /**< Bar container. */
   dictionary_type deathToBar_; /**< Map from death index to bar index. */
+  /**
+   * @brief Map from face ID to face position. Only stores a pair if ID != position.
+   */
+  std::unordered_map<id_index, pos_index> idToPosition_;  //TODO: test other map types
   bool isReduced_;              /**< True if `_reduce()` was called. */
 
   void _reduce();
@@ -114,19 +119,7 @@ class Base_pairing
 };
 
 template <class Master_matrix>
-inline Base_pairing<Master_matrix>::Base_pairing() : isReduced_(false) 
-{}
-
-template <class Master_matrix>
-inline Base_pairing<Master_matrix>::Base_pairing(const Base_pairing& matrixToCopy)
-    : barcode_(matrixToCopy.barcode_), deathToBar_(matrixToCopy.deathToBar_), isReduced_(matrixToCopy.isReduced_) 
-{}
-
-template <class Master_matrix>
-inline Base_pairing<Master_matrix>::Base_pairing(Base_pairing<Master_matrix>&& other) noexcept
-    : barcode_(std::move(other.barcode_)),
-      deathToBar_(std::move(other.deathToBar_)),
-      isReduced_(std::move(other.isReduced_)) 
+inline Base_pairing<Master_matrix>::Base_pairing() : PIDM(), isReduced_(false) 
 {}
 
 template <class Master_matrix>
@@ -139,7 +132,6 @@ inline const typename Base_pairing<Master_matrix>::barcode_type& Base_pairing<Ma
 template <class Master_matrix>
 inline void Base_pairing<Master_matrix>::_reduce() 
 {
-  using id_index = typename Master_matrix::index;
   std::unordered_map<id_index, index> pivotsToColumn(_matrix()->get_number_of_columns());
 
   auto dim = _matrix()->get_max_dimension();
@@ -175,8 +167,10 @@ inline void Base_pairing<Master_matrix>::_reduce()
 
         if (pivot != static_cast<id_index>(-1)) {
           pivotsToColumn.emplace(pivot, i);
-          _matrix()->get_column(pivot).clear();
-          barcode_.emplace_back(pivot, i, dim - 1);
+          auto it = idToPosition_.find(pivot);
+          auto pivotColumnNumber = it == idToPosition_.end() ? pivot : it->second;
+          _matrix()->get_column(pivotColumnNumber).clear();
+          barcode_.emplace_back(pivotColumnNumber, i, dim - 1);
         } else {
           curr.clear();
           barcode_.emplace_back(i, -1, dim);
@@ -216,15 +210,12 @@ inline void Base_pairing<Master_matrix>::_remove_last(pos_index columnIndex)
       deathToBar_.erase(it);
     };
   }
-}
 
-template <class Master_matrix>
-inline Base_pairing<Master_matrix>& Base_pairing<Master_matrix>::operator=(Base_pairing<Master_matrix> other) 
-{
-  barcode_.swap(other.barcode_);
-  deathToBar_.swap(other.deathToBar_);
-  std::swap(isReduced_, other.isReduced_);
-  return *this;
+  auto it = PIDM::map_.find(columnIndex);
+  if (it != PIDM::map_.end()){
+    idToPosition_.erase(it->second);
+    PIDM::map_.erase(it);
+  }
 }
 
 }  // namespace persistence_matrix

src/Persistence_matrix/include/gudhi/Persistence_matrix/boundary_face_position_to_id_mapper.h

@@ -0,0 +1,60 @@
+/*    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
+ */
+
+/**
+ * @file base_pairing.h
+ * @author Hannah Schreiber
+ * @brief Contains the @ref Gudhi::persistence_matrix::Base_pairing class and
+ * @ref Gudhi::persistence_matrix::Dummy_base_pairing structure.
+ */
+
+#ifndef PM_ID_POS_MAPPER_H
+#define PM_ID_POS_MAPPER_H
+
+#include <unordered_map>
+
+namespace Gudhi {
+namespace persistence_matrix {
+
+/**
+ * @private
+ * @ingroup persistence_matrix
+ *
+ * @brief Empty structure.
+ * Inherited instead of @ref Face_position_to_ID_mapper.
+ */
+struct Dummy_pos_mapper {
+  friend void swap([[maybe_unused]] Dummy_pos_mapper& d1, [[maybe_unused]] Dummy_pos_mapper& d2) {}
+};
+
+/**
+ * @private
+ * @ingroup persistence_matrix
+ *
+ * @brief Map from face position to face ID. Only stores a pair if ID != position and has_removable_column is true.
+ * 
+ * @tparam id_index @ref IDIdx index type
+ * @tparam pos_index @ref PosIdx index type
+ */
+template<typename id_index, typename pos_index>
+struct Face_position_to_ID_mapper {
+  using map_type = std::unordered_map<pos_index,id_index>; //TODO: test other map types
+
+  map_type map_;
+
+  friend void swap(Face_position_to_ID_mapper& mapper1, Face_position_to_ID_mapper& mapper2) {
+    mapper1.map_.swap(mapper2.map_);
+  }
+};
+
+}  // namespace persistence_matrix
+}  // namespace Gudhi
+
+#endif  // PM_ID_POS_MAPPER_H

src/Persistence_matrix/include/gudhi/Persistence_matrix/boundary_matrix.h

@@ -474,14 +474,16 @@ inline typename Boundary_matrix<Master_matrix>::index Boundary_matrix<Master_mat
 
   //updates container sizes
   if constexpr (Master_matrix::Option_list::has_row_access && !Master_matrix::Option_list::has_removable_rows) {
-    id_index pivot;
-    if constexpr (Master_matrix::Option_list::is_z2) {
-      pivot = *std::prev(boundary.end());
-    } else {
-      pivot = std::prev(boundary.end())->first;
+    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 container
-    if (ra_opt::rows_->size() <= pivot) ra_opt::rows_->resize(pivot + 1);
   }
 
   //row swap map containers
@@ -499,6 +501,16 @@ inline typename Boundary_matrix<Master_matrix>::index Boundary_matrix<Master_mat
     }
   }
 
+  //maps for possible shifting between column content and position indices used for birth events
+  if constexpr (activePairingOption){
+    if (faceIndex != nextInsertIndex_){
+      pair_opt::idToPosition_.emplace(faceIndex, nextInsertIndex_);
+      if constexpr (Master_matrix::Option_list::has_removable_columns){
+        pair_opt::PIDM::map_.emplace(nextInsertIndex_, faceIndex);
+      }
+    }
+  }
+
   _container_insert(boundary, nextInsertIndex_, dim);
 
   return nextInsertIndex_++;

src/Persistence_matrix/include/gudhi/Persistence_matrix/chain_vine_swap.h

@@ -84,6 +84,7 @@ class Chain_barcode_swap : public Chain_pairing<Master_matrix>
  public:
   using id_index = typename Master_matrix::id_index;    /**< @ref IDIdx index type. */
   using pos_index = typename Master_matrix::pos_index;  /**< @ref PosIdx index type. */
+  //CP = Chain Pairing
   using CP = Chain_pairing<Master_matrix>;
 
   /**

src/Persistence_matrix/include/gudhi/Persistence_matrix/overlay_ididx_to_matidx.h

@@ -622,6 +622,7 @@ class Id_to_index_overlay
 
   void _initialize_map(unsigned int size);
   index _id_to_index(id_index id) const;
+  index& _id_to_index(id_index id);
 };
 
 template <class Matrix_type, class Master_matrix_type>
@@ -640,7 +641,7 @@ inline Id_to_index_overlay<Matrix_type, Master_matrix_type>::Id_to_index_overlay
   _initialize_map(orderedBoundaries.size());
   if constexpr (Master_matrix_type::Option_list::is_of_boundary_type) {
     for (unsigned int i = 0; i < orderedBoundaries.size(); i++) {
-      idToIndex_->operator[](i) = i;
+      _id_to_index(i) = i;
     }
   }
 }
@@ -678,7 +679,7 @@ inline Id_to_index_overlay<Matrix_type, Master_matrix_type>::Id_to_index_overlay
   _initialize_map(orderedBoundaries.size());
   if constexpr (Master_matrix_type::Option_list::is_of_boundary_type) {
     for (unsigned int i = 0; i < orderedBoundaries.size(); i++) {
-      idToIndex_->operator[](i) = i;
+      _id_to_index(i) = i;
     }
   }
 }
@@ -739,7 +740,7 @@ inline void Id_to_index_overlay<Matrix_type, Master_matrix_type>::insert_boundar
       if (idToIndex_->size() == nextIndex_) {
         idToIndex_->push_back(nextIndex_);
       } else {
-        idToIndex_->operator[](nextIndex_) = nextIndex_;
+        _id_to_index(nextIndex_) = nextIndex_;
       }
     }
     ++nextIndex_;
@@ -756,7 +757,7 @@ inline void Id_to_index_overlay<Matrix_type, Master_matrix_type>::insert_boundar
     GUDHI_CHECK(idToIndex_->find(faceIndex) == idToIndex_->end(),
                 std::invalid_argument("Id_to_index_overlay::insert_boundary - Index for simplex already chosen!"));
   } else {
-    GUDHI_CHECK((idToIndex_->size() <= faceIndex || idToIndex_[faceIndex] == static_cast<index>(-1)),
+    GUDHI_CHECK((idToIndex_->size() <= faceIndex || _id_to_index(faceIndex) == static_cast<index>(-1)),
                 std::invalid_argument("Id_to_index_overlay::insert_boundary - Index for simplex already chosen!"));
   }
   matrix_.insert_boundary(faceIndex, boundary, dim);
@@ -767,7 +768,7 @@ inline void Id_to_index_overlay<Matrix_type, Master_matrix_type>::insert_boundar
       if (idToIndex_->size() <= faceIndex) {
         idToIndex_->resize(faceIndex + 1, -1);
       }
-      idToIndex_->operator[](faceIndex) = nextIndex_;
+      _id_to_index(faceIndex) = nextIndex_;
     }
     ++nextIndex_;
   }
@@ -804,7 +805,7 @@ inline void Id_to_index_overlay<Matrix_type, Master_matrix_type>::remove_maximal
       }
     } else {
       for (id_index i = 0; i < idToIndex_->size(); ++i) {
-        if (idToIndex_->operator[](i) != static_cast<index>(-1)) indexToID[idToIndex_->operator[](i)] = i;
+        if (_id_to_index(i) != static_cast<index>(-1)) indexToID[_id_to_index(i)] = i;
       }
     }
     --nextIndex_;
@@ -819,7 +820,7 @@ inline void Id_to_index_overlay<Matrix_type, Master_matrix_type>::remove_maximal
     if constexpr (Master_matrix_type::Option_list::has_map_column_container) {
       idToIndex_->erase(faceID);
     } else {
-      idToIndex_->operator[](faceID) = -1;
+      _id_to_index(faceID) = -1;
     }
   } else {
     matrix_.remove_maximal_face(faceID);
@@ -853,10 +854,10 @@ inline void Id_to_index_overlay<Matrix_type, Master_matrix_type>::remove_last()
       idToIndex_->erase(it);
     } else {
       index id = idToIndex_->size() - 1;
-      while (idToIndex_->operator[](id) == static_cast<index>(-1)) --id;  // should always stop before reaching -1
-      GUDHI_CHECK(idToIndex_->operator[](id) == nextIndex_,
+      while (_id_to_index(id) == static_cast<index>(-1)) --id;  // should always stop before reaching -1
+      GUDHI_CHECK(_id_to_index(id) == nextIndex_,
                   std::logic_error("Id_to_index_overlay::remove_last - Indexation problem."));
-      idToIndex_->operator[](id) = -1;
+      _id_to_index(id) = -1;
     }
   }
 }
@@ -1088,6 +1089,13 @@ Id_to_index_overlay<Matrix_type, Master_matrix_type>::_id_to_index(id_index id)
   }
 }
 
+template <class Matrix_type, class Master_matrix_type>
+inline typename Id_to_index_overlay<Matrix_type, Master_matrix_type>::index&
+Id_to_index_overlay<Matrix_type, Master_matrix_type>::_id_to_index(id_index id)
+{
+  return idToIndex_->operator[](id);  //for maps, the entry is created if not existing as needed in the constructors
+}
+
 }  // namespace persistence_matrix
 }  // namespace Gudhi