zigzag_persistence.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):       Clément Maria
 *
 *    Copyright (C) 2021 Inria
 *
 *    Modification(s):
 *      - 2023/05 Hannah Schreiber: Rework of the interface, reorganization and debug
 *      - 2023/05 Hannah Schreiber: Addition of infinite bars
 *      - 2024/06 Hannah Schreiber: Separation of the zigzag algorithm from the filtration value management
 *      - YYYY/MM Author: Description of the modification
 */
 
#ifndef ZIGZAG_PERSISTENCE_H_
#define ZIGZAG_PERSISTENCE_H_
 
#include <cmath>
#include <set>
#include <utility>
#include <vector>
 
#include <boost/iterator/indirect_iterator.hpp>
#if BOOST_VERSION >= 108100
#include <boost/unordered/unordered_flat_map.hpp>   //don't exist for lower versions of boost
// #include <boost/unordered/unordered_map.hpp>
#else
#include <unordered_map>
#endif
 
#include <gudhi/Matrix.h>
 
namespace Gudhi {
namespace zigzag_persistence {
 
template <Gudhi::persistence_matrix::Column_types column_type>
struct Zigzag_matrix_options : Gudhi::persistence_matrix::Default_options<column_type, true> {
  static const bool has_row_access = true;
  static const bool has_column_pairings = false;
  static const bool has_vine_update = true;
  static const bool is_of_boundary_type = false;
  static const bool has_map_column_container = true;
  static const bool has_removable_columns = true;
  static const bool has_removable_rows = true;
};
 
struct Default_zigzag_options {
  using Internal_key = int;   
  using Dimension = int;      
  static const Gudhi::persistence_matrix::Column_types column_type =
      Gudhi::persistence_matrix::Column_types::NAIVE_VECTOR;
};
 
// TODO: add the possibility of something else than Z2. Which means that the possibility of vineyards without Z2
// also needs to be implemented. The theory needs to be done first.
template <class ZigzagOptions = Default_zigzag_options>
class Zigzag_persistence
{
 public:
  using Options = ZigzagOptions;                    
  using Index = typename Options::Internal_key;     
  using Dimension = typename Options::Dimension;    
 private:
#if BOOST_VERSION >= 108100
  using Birth_dictionary = boost::unordered_flat_map<Index, Index>;   
  // using Birth_dictionary = boost::unordered_map<Index, Index>;        /**< Dictionary type. */
#else
  using Birth_dictionary = std::unordered_map<Index, Index>;             
#endif
  using Matrix_options = Zigzag_matrix_options<Options::column_type>; 
  using Matrix = Gudhi::persistence_matrix::Matrix<Matrix_options>;   
  using Matrix_index = typename Matrix::Index;                        
  struct Birth_ordering {
    Birth_ordering() : birthToPos_(), maxBirthPos_(0), minBirthPos_(-1) {}
 
    void add_birth_forward(Index arrowNumber) {  // amortized constant time
      birthToPos_.emplace_hint(birthToPos_.end(), arrowNumber, maxBirthPos_);
      ++maxBirthPos_;
    }
    void add_birth_backward(Index arrowNumber) {  // amortized constant time
      birthToPos_.emplace_hint(birthToPos_.end(), arrowNumber, minBirthPos_);
      --minBirthPos_;
    }
 
    void remove_birth(Index birth) { birthToPos_.erase(birth); }
    bool birth_order(Index k1, Index k2) const { return birthToPos_.at(k1) < birthToPos_.at(k2); }
    bool reverse_birth_order(Index k1, Index k2) const { return birthToPos_.at(k1) > birthToPos_.at(k2); }
 
   private:
    Birth_dictionary birthToPos_; 
    Index maxBirthPos_;            
    Index minBirthPos_;            
  };
 
 public:
  Zigzag_persistence(std::function<void(Dimension, Index, Index)> stream_interval,
                     unsigned int preallocationSize = 0)
      : matrix_(
            preallocationSize,
            [this](Matrix_index columnIndex1, Matrix_index columnIndex2) -> bool {
              if (matrix_.get_column(columnIndex1).is_paired()) {
                return matrix_.get_pivot(columnIndex1) < matrix_.get_pivot(columnIndex2);
              }
              return birthOrdering_.birth_order(births_.at(columnIndex1), births_.at(columnIndex2));
            },
            [this](Matrix_index columnIndex1, Matrix_index columnIndex2) -> bool { return false; }),
        numArrow_(-1),
        stream_interval_(std::move(stream_interval)) {}
  template <class BoundaryRange = std::initializer_list<Index> >
  Index insert_cell(const BoundaryRange& boundary, Dimension dimension) {
    ++numArrow_;
    _process_forward_arrow(boundary, dimension);
    return numArrow_;
  }
 
  Index remove_cell(Index arrowNumber) {
    ++numArrow_;
    _process_backward_arrow(arrowNumber);
    return numArrow_;
  }
 
  Index apply_identity() { return ++numArrow_; }
 
  template <typename F>
  void get_current_infinite_intervals(F&& stream_infinite_interval) {
    for (auto& p : births_) {
      auto& col = matrix_.get_column(p.first);
      stream_infinite_interval(col.get_dimension(), p.second);
    }
  }
 
 private:
  template <class BoundaryRange>
  void _process_forward_arrow(const BoundaryRange& boundary, Dimension dim) {
    std::vector<Matrix_index> chainsInF = matrix_.insert_boundary(numArrow_, boundary, dim);
 
    if (!chainsInF.empty()) {
      _apply_surjective_reflection_diamond(dim, chainsInF);
    } else {
      birthOrdering_.add_birth_forward(numArrow_);
      births_[matrix_.get_column_with_pivot(numArrow_)] = numArrow_;
    }
  }
 
  void _apply_surjective_reflection_diamond(Dimension dim, const std::vector<Matrix_index>& chainsInF) {
    // fp is the largest death index for <=d
    // Set col_fp: col_fp <- col_f1+...+col_fp (now in G); preserves lowest idx
    auto chainFp = chainsInF[0];  // col_fp, with largest death <d index.
 
    // chainsInF is ordered, from .begin() to end(), by decreasing lowest_idx_. The
    // lowest_idx_ is also the death of the chain in the right suffix of the
    // filtration (all backward arrows). Consequently, the chains in F are ordered by
    // decreasing death for <d.
    // Pair the col_fi, i = 1 ... p-1, according to the reflection diamond principle
    // Order the fi by reverse birth ordering <=_b
    auto cmp_birth = [this](Index k1, Index k2) -> bool {
      return birthOrdering_.reverse_birth_order(k1, k2);
    };  // true iff b(k1) >b b(k2)
 
    // availableBirth: for all i by >d value of the d_i,
    // contains at step i all b_j, j > i, and maybe b_i if not stolen
    std::set<Index, decltype(cmp_birth)> availableBirth(cmp_birth);
    // for f1 to f_{p} (i by <=d), insertion in availableBirth sorts by >=b
    for (auto& chainF : chainsInF) {
      availableBirth.insert(births_.at(chainF));
    }
 
    auto maxbIt = availableBirth.begin();  // max birth cycle
    auto maxb = *maxbIt;                   // max birth value, for persistence diagram
    availableBirth.erase(maxbIt);          // remove max birth cycle (stolen)
 
    auto lastModifiedChainIt = chainsInF.rbegin();
 
    // consider all death indices by increasing <d order i.e., increasing lowest_idx_
    for (auto chainFIt = chainsInF.rbegin();  // by increasing death order <d
         *chainFIt != chainFp; ++chainFIt)    // chain_fp=*begin() has max death
    {                                         // find which reduced col has this birth
      auto birthIt = availableBirth.find(births_.at(*chainFIt));
      if (birthIt == availableBirth.end())  // birth is not available. *chain_f_it
      {                                     // must become the sum of all chains in F with smaller death index.
        // this gives as birth the maximal birth of all chains with strictly larger
        // death <=> the maximal available death.
        // Let c_1 ... c_f be the chains s.t. <[c_1+...+c_f]> is the kernel and
        //  death(c_i) >d death(c_i-1). If the birth of c_i is not available, we set
        // c_i <- c_i + c_i-1 + ... + c_1, which is [c_i + c_i-1 + ... + c_1] on
        // the right (of death the maximal<d death(c_i)), and is [c_i + c_i-1 + ... +
        // c_1] + kernel = [c_f + c_f-1 + ... + c_i+1] on the left (of birth the max<b
        // of the birth of the c_j, j>i  <=> the max<b available birth).
        // N.B. some of the c_k, k<i, have already been modified to be equal to
        // c_k + c_k-1 + ... + c_1. The largest k with this property is maintained in
        // last_modified_chain_it (no need to compute from scratch the full sum).
 
        // last_modified is equal to c_k+...+c_1, all c_j, i>j>k, are indeed c_j
        // set c_i <- c_i + (c_i-1) + ... + (c_k+1) + (c_k + ... + c_1)
        for (auto chainPassedIt = lastModifiedChainIt; chainPassedIt != chainFIt; ++chainPassedIt) {
          // all with smaller <d death
          matrix_.add_to(*chainPassedIt, *chainFIt);
        }
        lastModifiedChainIt = chainFIt;  // new cumulated c_i+...+c_1
        // remove the max available death
        auto maxAvailBirthIt = availableBirth.begin();  // max because order by decreasing <b
        Index maxAvailBirth = *maxAvailBirthIt;         // max available birth
 
        births_.at(*chainFIt) = maxAvailBirth;  // give new birth
        availableBirth.erase(maxAvailBirthIt);  // remove birth from availability
      } else {
        availableBirth.erase(birthIt);
      }  // birth not available anymore, do not
    }  // modify *chain_f_it.
 
    birthOrdering_.remove_birth(maxb);
    births_.erase(chainFp);
 
    // Update persistence diagram with left interval [fil(b_max) ; fil(m))
    stream_interval_(dim - 1, maxb, numArrow_);
  }
 
  void _process_backward_arrow(Index cellID) {
    // column whose key is the one of the removed cell
    Matrix_index currCol = matrix_.get_column_with_pivot(cellID);
 
    // Record all columns that get affected by the transpositions, i.e., have a coeff
    std::vector<Matrix_index> modifiedColumns;
    const auto& row = matrix_.get_row(cellID);
    modifiedColumns.reserve(row.size());
    std::transform(row.begin(), row.end(), std::back_inserter(modifiedColumns),
                   [](const auto& cell) { return cell.get_column_index(); });
    // Sort by left-to-right order in the matrix_ (no order maintained in rows)
    std::stable_sort(modifiedColumns.begin(), modifiedColumns.end(), [this](Matrix_index i1, Matrix_index i2) {
      return matrix_.get_pivot(i1) < matrix_.get_pivot(i2);
    });
 
    // Modifies curr_col, not the other one.
    for (auto otherColIt = std::next(modifiedColumns.begin()); otherColIt != modifiedColumns.end(); ++otherColIt) {
      currCol = matrix_.vine_swap_with_z_eq_1_case(currCol, *otherColIt);
    }
 
    // curr_col points to the column to remove by restriction of K to K-{\sigma}
    auto& col = matrix_.get_column(currCol);
    if (!col.is_paired()) {  // in F
      auto it = births_.find(currCol);
      stream_interval_(col.get_dimension(), it->second, numArrow_);
      birthOrdering_.remove_birth(it->second);
      births_.erase(it);
    } else {  // in H    -> paired with c_g, that now belongs to F now
      // maintain the <=b order
      birthOrdering_.add_birth_backward(numArrow_);
      births_[col.get_paired_chain_index()] = numArrow_;
    }
 
    // cannot be in G as the removed cell is maximal
    matrix_.remove_maximal_cell(cellID, {});  // also un-pairs c_g if in H
  }
 
 private:
  Matrix matrix_;           
  Birth_dictionary births_;      
  Birth_ordering birthOrdering_; 
  Index numArrow_;               
  std::function<void(Dimension, Index, Index)> stream_interval_; 
};  // end class Zigzag_persistence
 
}  // namespace zigzag_persistence
}  // namespace Gudhi
 
#endif  // ZIGZAG_PERSISTENCE_H_