Loading...
Searching...
No Matches
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>
#include <boost/version.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_
Matrix.h
Contains Gudhi::persistence_matrix::Matrix class.
Gudhi::persistence_matrix::Matrix
Data structure for matrices, and in particular thought for matrices representing filtered complexes i...
Definition Matrix.h:144
Gudhi::persistence_matrix::Matrix< Matrix_options >::Index
typename Matrix_options::Index Index
Definition Matrix.h:147
Gudhi::zigzag_persistence::Zigzag_persistence::Zigzag_persistence
Zigzag_persistence(std::function< void(Dimension, Index, Index)> stream_interval, unsigned int preallocationSize=0)
Constructor of the Zigzag_persistence class.
Definition zigzag_persistence.h:264
Gudhi::zigzag_persistence::Zigzag_persistence::remove_cell
Index remove_cell(Index arrowNumber)
Updates the zigzag persistence diagram after the removal of the given cell.
Definition zigzag_persistence.h:301
Gudhi::zigzag_persistence::Zigzag_persistence::insert_cell
Index insert_cell(const BoundaryRange &boundary, Dimension dimension)
Updates the zigzag persistence diagram after the insertion of the given cell.
Definition zigzag_persistence.h:289
Gudhi::zigzag_persistence::Zigzag_persistence::apply_identity
Index apply_identity()
To use when a cell is neither inserted nor removed, but the filtration moves along the identity opera...
Definition zigzag_persistence.h:313
Gudhi::zigzag_persistence::Zigzag_persistence::Dimension
typename Options::Dimension Dimension
Definition zigzag_persistence.h:154
Gudhi::zigzag_persistence::Zigzag_persistence::Options
ZigzagOptions Options
Definition zigzag_persistence.h:152
Gudhi::zigzag_persistence::Zigzag_persistence::Index
typename Options::Internal_key Index
Definition zigzag_persistence.h:153
Gudhi::zigzag_persistence::Zigzag_persistence::get_current_infinite_intervals
void get_current_infinite_intervals(F &&stream_infinite_interval)
Outputs through the given callback method all birth indices which are currently not paired with a dea...
Definition zigzag_persistence.h:324
Gudhi::persistence_matrix::Column_types
Column_types
List of column types.
Definition persistence_matrix_options.h:30
Gudhi::persistence_matrix::Column_types::NAIVE_VECTOR
@ NAIVE_VECTOR
Definition persistence_matrix_options.h:37
Gudhi
Gudhi namespace.
Definition SimplicialComplexForAlpha.h:14
Gudhi::persistence_matrix::Default_options
Default option structure for Matrix class. See the PersistenceMatrixOptions concept for a more detail...
Definition persistence_matrix_options.h:79
Gudhi::zigzag_persistence::Default_zigzag_options
Default options for Zigzag_persistence.
Definition zigzag_persistence.h:65
Gudhi::zigzag_persistence::Default_zigzag_options::Internal_key
int Internal_key
Definition zigzag_persistence.h:66
Gudhi::zigzag_persistence::Default_zigzag_options::Dimension
int Dimension
Definition zigzag_persistence.h:67
Gudhi::zigzag_persistence::Default_zigzag_options::column_type
static const Gudhi::persistence_matrix::Column_types column_type
Column type use by the internal matrix.
Definition zigzag_persistence.h:71
Gudhi::zigzag_persistence::Zigzag_matrix_options
Options for the internal matrix of Zigzag_persistence.
Definition zigzag_persistence.h:50
Gudhi::zigzag_persistence::ZigzagOptions
List of options used for the zigzag persistence computation.
Definition ZigzagOptions.h:60
Gudhi::zigzag_persistence::ZigzagOptions::Dimension
unspecified Dimension
Type for the dimension values.
Definition ZigzagOptions.h:69
Gudhi::zigzag_persistence::ZigzagOptions::Internal_key
unspecified Internal_key
Numerical type for the cell IDs used internally and other indexations. It must be signed.
Definition ZigzagOptions.h:64