filtered_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 FILTERED_ZIGZAG_PERSISTENCE_H_
#define FILTERED_ZIGZAG_PERSISTENCE_H_
 
#include <cmath>
#include <limits>
#include <unordered_map>
#include <utility>
#include <vector>
 
#include <gudhi/Debug_utils.h>
#include <gudhi/zigzag_persistence.h>
#include <gudhi/persistence_interval.h>
 
namespace Gudhi {
namespace zigzag_persistence {
 
struct Default_filtered_zigzag_options : Default_zigzag_options {
  using Cell_key = int;            
  using Filtration_value = double; 
};
 
template <class FilteredZigzagOptions = Default_filtered_zigzag_options>
class Filtered_zigzag_persistence_with_storage
{
 public:
  using Options = FilteredZigzagOptions;                        
  using Internal_key = typename Options::Internal_key;          
  using Cell_key = typename Options::Cell_key;                  
  using Filtration_value = typename Options::Filtration_value;  
  using Dimension = typename Options::Dimension;                
  using Index_interval = Gudhi::persistence_matrix::Persistence_interval<Dimension,Internal_key>;
  using Filtration_value_interval = Gudhi::persistence_matrix::Persistence_interval<Dimension,Filtration_value>;
 
  Filtered_zigzag_persistence_with_storage(unsigned int preallocationSize = 0, int ignoreCyclesAboveDim = -1)
      : dimMax_(ignoreCyclesAboveDim),
        persistenceDiagram_(),
        numArrow_(-1),
        previousFiltrationValue_(std::numeric_limits<Filtration_value>::infinity()),
        pers_(
            [&](Dimension dim, Internal_key birth, Internal_key death) {
              if (dimMax_ == -1 || (dimMax_ != -1 && dim < dimMax_)) {  // don't record intervals over max dim
                persistenceDiagram_.emplace_back(birth, death, dim);
              }
            },
            preallocationSize) {}
 
  template <class BoundaryRange = std::initializer_list<Cell_key> >
  Internal_key insert_cell(Cell_key cellID,
                           const BoundaryRange& boundary,
                           Dimension dimension,
                           Filtration_value filtrationValue)
  {
    if (dimMax_ != -1 && dimension > dimMax_) {
      return apply_identity();
    }
 
    ++numArrow_;
 
    _store_filtration_value(filtrationValue);
 
    [[maybe_unused]] auto res = handleToKey_.try_emplace(cellID, numArrow_);
 
    GUDHI_CHECK(res.second, "Zigzag_persistence::insert_cell - cell already in the complex");
 
    // Compute the keys of the cells of the boundary.
    std::vector<Internal_key> translatedBoundary;
    translatedBoundary.reserve(dimension * 2);  // boundary does not have to have `size()`
    for (auto b : boundary) {
      translatedBoundary.push_back(handleToKey_.at(b));  // TODO: add possibilities of coefficients
    }
    std::sort(translatedBoundary.begin(), translatedBoundary.end());
 
    pers_.insert_cell(translatedBoundary, dimension);
 
    return numArrow_;
  }
 
  Internal_key remove_cell(Cell_key cellID, Filtration_value filtrationValue) {
    auto it = handleToKey_.find(cellID);
 
    if (it == handleToKey_.end()) {
      return apply_identity();
    }
 
    ++numArrow_;
 
    _store_filtration_value(filtrationValue);
 
    pers_.remove_cell(it->second);
    handleToKey_.erase(it);
 
    return numArrow_;
  }
 
  Internal_key apply_identity() {
    ++numArrow_;
    pers_.apply_identity();
    return numArrow_;
  }
 
  const std::vector<Index_interval>& get_index_persistence_diagram() const { return persistenceDiagram_; }
 
  Filtration_value get_filtration_value_from_index(Internal_key idx) {
    // lower_bound(x) returns leftmost y s.t. x <= y
    auto itBirth =
        std::lower_bound(filtrationValues_.begin(),
                         filtrationValues_.end(),
                         idx,
                         [](std::pair<Internal_key, Filtration_value> p, Internal_key k) { return p.first < k; });
    if (itBirth == filtrationValues_.end() || itBirth->first > idx) {
      --itBirth;
    }
    return itBirth->second;
  };
 
  std::vector<Filtration_value_interval> get_persistence_diagram(Filtration_value shortestInterval = 0.,
                                                                 bool includeInfiniteBars = true) {
    std::vector<Filtration_value_interval> diag = _get_persistence_diagram(shortestInterval);
 
    if (includeInfiniteBars) {
      _retrieve_infinite_bars(diag);
    }
 
    return diag;
  }
 
 private:
  std::unordered_map<Cell_key, Internal_key> handleToKey_;  
  Dimension dimMax_;                                        
  std::vector<Index_interval> persistenceDiagram_;          
  Internal_key numArrow_;                                   
  Filtration_value previousFiltrationValue_;                
  std::vector<std::pair<Internal_key, Filtration_value> > filtrationValues_;
  Zigzag_persistence<FilteredZigzagOptions> pers_;          
  void _store_filtration_value(Filtration_value filtrationValue) {
    if (filtrationValue != previousFiltrationValue_)  // check whether the filt value has changed
    {
      // consecutive pairs (i,f), (j,f') mean cells of index k in [i,j-1] have filtration value f
      previousFiltrationValue_ = filtrationValue;
      filtrationValues_.emplace_back(numArrow_, previousFiltrationValue_);
    }
  }
 
  std::vector<Filtration_value_interval> _get_persistence_diagram(Filtration_value shortestInterval) {
    std::vector<Filtration_value_interval> diag;
    diag.reserve(persistenceDiagram_.size());
 
    // std::stable_sort(filtrationValues_.begin(), filtrationValues_.end(),
    //                  [](std::pair<Internal_key, Filtration_value> p1, std::pair<Internal_key, Filtration_value> p2) {
    //                    return p1.first < p2.first;
    //                  });
 
    for (auto bar : persistenceDiagram_) {
      Filtration_value birth = get_filtration_value_from_index(bar.birth);
      Filtration_value death = get_filtration_value_from_index(bar.death);
      if (birth > death) {
        std::swap(birth, death);
      }
 
      if (death - birth > shortestInterval) {
        diag.emplace_back(birth, death, bar.dim);
      }
    }
 
    return diag;
  }
 
  void _retrieve_infinite_bars(std::vector<Filtration_value_interval>& diag) {
    auto stream_infinite_interval = [&](Dimension dim, Internal_key birthIndex) {
      if (dimMax_ == -1 || (dimMax_ != -1 && dim < dimMax_))
        diag.emplace_back(get_filtration_value_from_index(birthIndex), Filtration_value_interval::inf, dim);
    };
 
    pers_.get_current_infinite_intervals(stream_infinite_interval);
  }
};  // end class Filtered_zigzag_persistence_with_storage
 
template <class FilteredZigzagOptions = Default_filtered_zigzag_options>
class Filtered_zigzag_persistence {
 public:
  using Options = FilteredZigzagOptions;                        
  using Internal_key = typename Options::Internal_key;          
  using Cell_key = typename Options::Cell_key;                  
  using Filtration_value = typename Options::Filtration_value;  
  using Dimension = typename Options::Dimension;      
  template <typename F>
  Filtered_zigzag_persistence(F&& stream_interval, unsigned int preallocationSize = 0)
      : handleToKey_(preallocationSize),
        numArrow_(-1),
        keyToFiltrationValue_(preallocationSize),
        pers_(
            [&, stream_interval = std::forward<F>(stream_interval)](Dimension dim,
                                                                    Internal_key birth,
                                                                    Internal_key death) {
              auto itB = keyToFiltrationValue_.find(birth);
              auto itD = keyToFiltrationValue_.find(death);
              if (itB->second != itD->second) stream_interval(dim, itB->second, itD->second);
              keyToFiltrationValue_.erase(itB);
              keyToFiltrationValue_.erase(itD);
            },
            preallocationSize) {}
 
  template <class BoundaryRange = std::initializer_list<Cell_key> >
  Internal_key insert_cell(Cell_key cellID,
                           const BoundaryRange& boundary,
                           Dimension dimension,
                           Filtration_value filtrationValue)
  {
    ++numArrow_;
 
    [[maybe_unused]] auto res = handleToKey_.try_emplace(cellID, numArrow_);
 
    GUDHI_CHECK(res.second, "Zigzag_persistence::insert_cell - cell already in the complex");
 
    keyToFiltrationValue_.try_emplace(numArrow_, filtrationValue);
 
    // Compute the keys of the cells of the boundary.
    std::vector<Internal_key> translatedBoundary;
    translatedBoundary.reserve(dimension * 2);  // boundary does not have to have `size()`
    for (auto b : boundary) {
      translatedBoundary.push_back(handleToKey_.at(b));  // TODO: add possibilities of coefficients
    }
    std::sort(translatedBoundary.begin(), translatedBoundary.end());
 
    pers_.insert_cell(translatedBoundary, dimension);
 
    return numArrow_;
  }
 
  Internal_key remove_cell(Cell_key cellID, Filtration_value filtrationValue) {
    ++numArrow_;
 
    auto it = handleToKey_.find(cellID);
    GUDHI_CHECK(it != handleToKey_.end(), "Zigzag_persistence::remove_cell - cell not in the complex");
 
    keyToFiltrationValue_.try_emplace(numArrow_, filtrationValue);
 
    pers_.remove_cell(it->second);
    handleToKey_.erase(it);
 
    return numArrow_;
  }
 
  Internal_key apply_identity() {
    ++numArrow_;
    pers_.apply_identity();
    return numArrow_;
  }
 
  template <typename F>
  void get_current_infinite_intervals(F&& stream_infinite_interval) {
    pers_.get_current_infinite_intervals(
        [&](Dimension dim, Internal_key birth) { stream_infinite_interval(dim, keyToFiltrationValue_.at(birth)); });
  }
 
 private:
  template <typename key_type, typename value_type>
  using Dictionary = std::unordered_map<key_type, value_type>;  // TODO: benchmark with other map types
 
  Dictionary<Cell_key, Internal_key> handleToKey_;                  
  Internal_key numArrow_;                                           
  Dictionary<Internal_key, Filtration_value> keyToFiltrationValue_; 
  Zigzag_persistence<FilteredZigzagOptions> pers_;                  
};  // end class Filtered_zigzag_persistence
 
}  // namespace zigzag_persistence
}  // namespace Gudhi
 
#endif  // FILTERED_ZIGZAG_PERSISTENCE_H_