text: /*    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) 2014 Inria
   *
   *    Modification(s):
   *      - YYYY/MM Author: Description of the modification
   */
  
  #ifndef PERSISTENT_COHOMOLOGY_H_
  #define PERSISTENT_COHOMOLOGY_H_
  
  #include <gudhi/Persistent_cohomology/Persistent_cohomology_column.h>
  #include <gudhi/Persistent_cohomology/Field_Zp.h>
  #include <gudhi/Simple_object_pool.h>
  
  #include <boost/intrusive/set.hpp>
  #include <boost/pending/disjoint_sets.hpp>
  #include <boost/intrusive/list.hpp>
  
  #include <map>
  #include <utility>
  #include <list>
  #include <vector>
  #include <set>
  #include <fstream>  // std::ofstream
  #include <limits>  // for numeric_limits<>
  #include <tuple>
  #include <algorithm>
  #include <string>
  #include <stdexcept>  // for std::out_of_range
  
  namespace Gudhi {
  
  namespace persistent_cohomology {
  
  /** \brief Computes the persistent cohomology of a filtered complex.
   *
   * \ingroup persistent_cohomology
   * 
   * The computation is implemented with a Compressed Annotation Matrix
   * (CAM)\cite DBLP:conf/esa/BoissonnatDM13,
   * and is adapted to the computation of Multi-Field Persistent Homology (MF)
   * \cite boissonnat:hal-00922572 .
   *
   * \implements PersistentHomology
   *
   */
  // TODO(CM): Memory allocation policy: classic, use a mempool, etc.
  template<class FilteredComplex, class CoefficientField>
  class Persistent_cohomology {
   public:
    // Data attached to each simplex to interface with a Property Map.
  
    /** \brief Data stored for each simplex. */
    typedef typename FilteredComplex::Simplex_key Simplex_key;
    /** \brief Handle to specify a simplex. */
    typedef typename FilteredComplex::Simplex_handle Simplex_handle;
    /** \brief Type for the value of the filtration function. */
    typedef typename FilteredComplex::Filtration_value Filtration_value;
    /** \brief Type of element of the field. */
    typedef typename CoefficientField::Element Arith_element;
    /** \brief Type for birth and death FilteredComplex::Simplex_handle.
     * The Arith_element field is used for the multi-field framework. */
    typedef std::tuple<Simplex_handle, Simplex_handle, Arith_element> Persistent_interval;
  
   private:
    // Compressed Annotation Matrix types:
    // Column type
    typedef Persistent_cohomology_column<Simplex_key, Arith_element> Column;  // contains 1 set_hook
    // Cell type
    typedef typename Column::Cell Cell;   // contains 2 list_hooks
    // Remark: constant_time_size must be false because base_hook_cam_h has auto_unlink link_mode
    typedef boost::intrusive::list<Cell,
        boost::intrusive::constant_time_size<false>,
        boost::intrusive::base_hook<base_hook_cam_h> > Hcell;
  
    typedef boost::intrusive::set<Column,
        boost::intrusive::constant_time_size<false> > Cam;
    // Sparse column type for the annotation of the boundary of an element.
    typedef std::vector<std::pair<Simplex_key, Arith_element> > A_ds_type;
  
   public:
    /** \brief Initializes the Persistent_cohomology class.
     *
     * @param[in] cpx Complex for which the persistent homology is computed.
     * cpx is a model of FilteredComplex
     *
     * @param[in] persistence_dim_max if true, the persistent homology for the maximal dimension in the
     *                                complex is computed. If false, it is ignored. Default is false.
     *
     * @exception std::out_of_range In case the number of simplices is more than Simplex_key type numeric limit.
     */
    explicit Persistent_cohomology(FilteredComplex& cpx, bool persistence_dim_max = false)
        : cpx_(&cpx),
          dim_max_(cpx.dimension()),                       // upper bound on the dimension of the simplices
          coeff_field_(),                                  // initialize the field coefficient structure.
          num_simplices_(cpx_->num_simplices()),           // num_simplices save to avoid to call thrice the function
          ds_rank_(num_simplices_),                        // union-find
          ds_parent_(num_simplices_),                      // union-find
          ds_repr_(num_simplices_, NULL),                  // union-find -> annotation vectors
          dsets_(&ds_rank_[0], &ds_parent_[0]),            // union-find
          cam_(),                                          // collection of annotation vectors
          zero_cocycles_(),                                // union-find -> Simplex_key of creator for 0-homology
          transverse_idx_(),                               // key -> row
          persistent_pairs_(),
          interval_length_policy(&cpx, 0),
          column_pool_(),  // memory pools for the CAM
          cell_pool_() {
      if (cpx_->num_simplices() > std::numeric_limits<Simplex_key>::max()) {
        // num_simplices must be strictly lower than the limit, because a value is reserved for null_key.
        throw std::out_of_range("The number of simplices is more than Simplex_key type numeric limit.");
      }
      Simplex_key idx_fil = 0;
      for (auto sh : cpx_->filtration_simplex_range()) {
        cpx_->assign_key(sh, idx_fil);
        ++idx_fil;
        dsets_.make_set(cpx_->key(sh));
      }
      if (persistence_dim_max) {
        ++dim_max_;
      }
    }
  
    ~Persistent_cohomology() {
      // Clean the transversal lists
      for (auto & transverse_ref : transverse_idx_) {
        // Destruct all the cells
        transverse_ref.second.row_->clear_and_dispose([&](Cell*p){p->~Cell();});
        delete transverse_ref.second.row_;
      }
    }
  
   private:
    struct length_interval {
      length_interval(FilteredComplex * cpx, Filtration_value min_length)
          : cpx_(cpx),
            min_length_(min_length) {
      }
  
      bool operator()(Simplex_handle sh1, Simplex_handle sh2) {
        return cpx_->filtration(sh2) - cpx_->filtration(sh1) > min_length_;
      }
  
      void set_length(Filtration_value new_length) {
        min_length_ = new_length;
      }
  
      FilteredComplex * cpx_;
      Filtration_value min_length_;
    };
  
   public:
    /** \brief Initializes the coefficient field.*/
    void init_coefficients(int charac) {
      coeff_field_.init(charac);
    }
    /** \brief Initializes the coefficient field for multi-field persistent homology.*/
    void init_coefficients(int charac_min, int charac_max) {
      coeff_field_.init(charac_min, charac_max);
    }
  
    /** \brief Compute the persistent homology of the filtered simplicial
     * complex.
     *
     * @param[in] min_interval_length the computation discards all intervals of length
     *                                less or equal than min_interval_length
     *
     * Assumes that the filtration provided by the simplicial complex is
     * valid. Undefined behavior otherwise. */
    void compute_persistent_cohomology(Filtration_value min_interval_length = 0) {
      interval_length_policy.set_length(min_interval_length);
      // Compute all finite intervals
      for (auto sh : cpx_->filtration_simplex_range()) {
        int dim_simplex = cpx_->dimension(sh);
        switch (dim_simplex) {
          case 0:
            break;
          case 1:
            update_cohomology_groups_edge(sh);
            break;
          default:
            update_cohomology_groups(sh, dim_simplex);
            break;
        }
      }
      // Compute infinite intervals of dimension 0
      Simplex_key key;
      for (auto v_sh : cpx_->skeleton_simplex_range(0)) {  // for all 0-dimensional simplices
        key = cpx_->key(v_sh);
  
        if (ds_parent_[key] == key  // root of its tree
        && zero_cocycles_.find(key) == zero_cocycles_.end()) {
          persistent_pairs_.emplace_back(
              cpx_->simplex(key), cpx_->null_simplex(), coeff_field_.characteristic());
        }
      }
      for (auto zero_idx : zero_cocycles_) {
        persistent_pairs_.emplace_back(
            cpx_->simplex(zero_idx.second), cpx_->null_simplex(), coeff_field_.characteristic());
      }
      // Compute infinite interval of dimension > 0
      for (auto cocycle : transverse_idx_) {
        persistent_pairs_.emplace_back(
            cpx_->simplex(cocycle.first), cpx_->null_simplex(), cocycle.second.characteristics_);
      }
    }
  
   private:
    /** \brief Update the cohomology groups under the insertion of an edge.
     *
     * The 0-homology is maintained with a simple Union-Find data structure, which
     * explains the existance of a specific function of edge insertions. */
    void update_cohomology_groups_edge(Simplex_handle sigma) {
      Simplex_handle u, v;
      boost::tie(u, v) = cpx_->endpoints(sigma);
  
      Simplex_key ku = dsets_.find_set(cpx_->key(u));
      Simplex_key kv = dsets_.find_set(cpx_->key(v));
  
      if (ku != kv) {        // Destroy a connected component
        dsets_.link(ku, kv);
        // Keys of the simplices which created the connected components containing
        // respectively u and v.
        Simplex_key idx_coc_u, idx_coc_v;
        auto map_it_u = zero_cocycles_.find(ku);
        // If the index of the cocycle representing the class is already ku.
        if (map_it_u == zero_cocycles_.end()) {
          idx_coc_u = ku;
        } else {
          idx_coc_u = map_it_u->second;
        }
  
        auto map_it_v = zero_cocycles_.find(kv);
        // If the index of the cocycle representing the class is already kv.
        if (map_it_v == zero_cocycles_.end()) {
          idx_coc_v = kv;
        } else {
          idx_coc_v = map_it_v->second;
        }
  
        if (cpx_->filtration(cpx_->simplex(idx_coc_u))
            < cpx_->filtration(cpx_->simplex(idx_coc_v))) {  // Kill cocycle [idx_coc_v], which is younger.
          if (interval_length_policy(cpx_->simplex(idx_coc_v), sigma)) {
            persistent_pairs_.emplace_back(
                cpx_->simplex(idx_coc_v), sigma, coeff_field_.characteristic());
          }
          // Maintain the index of the 0-cocycle alive.
          if (kv != idx_coc_v) {
            zero_cocycles_.erase(map_it_v);
          }
          if (kv == dsets_.find_set(kv)) {
            if (ku != idx_coc_u) {
              zero_cocycles_.erase(map_it_u);
            }
            zero_cocycles_[kv] = idx_coc_u;
          }
        } else {  // Kill cocycle [idx_coc_u], which is younger.
          if (interval_length_policy(cpx_->simplex(idx_coc_u), sigma)) {
            persistent_pairs_.emplace_back(
                cpx_->simplex(idx_coc_u), sigma, coeff_field_.characteristic());
          }
          // Maintain the index of the 0-cocycle alive.
          if (ku != idx_coc_u) {
            zero_cocycles_.erase(map_it_u);
          }
          if (ku == dsets_.find_set(ku)) {
            if (kv != idx_coc_v) {
              zero_cocycles_.erase(map_it_v);
            }
            zero_cocycles_[ku] = idx_coc_v;
          }
        }
        cpx_->assign_key(sigma, cpx_->null_key());
      } else if (dim_max_ > 1) {  // If ku == kv, same connected component: create a 1-cocycle class.
        create_cocycle(sigma, coeff_field_.multiplicative_identity(), coeff_field_.characteristic());
      }
    }
  
    /*
     * Compute the annotation of the boundary of a simplex.
     */
    void annotation_of_the_boundary(
        std::map<Simplex_key, Arith_element> & map_a_ds, Simplex_handle sigma,
        int dim_sigma) {
      // traverses the boundary of sigma, keeps track of the annotation vectors,
      // with multiplicity. We used to sum the coefficients directly in
      // annotations_in_boundary by using a map, we now do it later.
      typedef std::pair<Column *, int> annotation_t;
  #ifdef GUDHI_CAN_USE_CXX11_THREAD_LOCAL
      thread_local
  #endif  // GUDHI_CAN_USE_CXX11_THREAD_LOCAL
      std::vector<annotation_t> annotations_in_boundary;
      annotations_in_boundary.clear();
      int sign = 1 - 2 * (dim_sigma % 2);  // \in {-1,1} provides the sign in the
                                           // alternate sum in the boundary.
      Simplex_key key;
      Column * curr_col;
  
      for (auto sh : cpx_->boundary_simplex_range(sigma)) {
        key = cpx_->key(sh);
        if (key != cpx_->null_key()) {  // A simplex with null_key is a killer, and have null annotation
          // Find its annotation vector
          curr_col = ds_repr_[dsets_.find_set(key)];
          if (curr_col != NULL) {  // and insert it in annotations_in_boundary with multyiplicative factor "sign".
            annotations_in_boundary.emplace_back(curr_col, sign);
          }
        }
        sign = -sign;
      }
      // Place identical annotations consecutively so we can easily sum their multiplicities.
      std::sort(annotations_in_boundary.begin(), annotations_in_boundary.end(),
                [](annotation_t const& a, annotation_t const& b) { return a.first < b.first; });
  
      // Sum the annotations with multiplicity, using a map<key,coeff>
      // to represent a sparse vector.
      std::pair<typename std::map<Simplex_key, Arith_element>::iterator, bool> result_insert_a_ds;
  
      for (auto ann_it = annotations_in_boundary.begin(); ann_it != annotations_in_boundary.end(); /**/) {
        Column* col = ann_it->first;
        int mult = ann_it->second;
        while (++ann_it != annotations_in_boundary.end() && ann_it->first == col) {
          mult += ann_it->second;
        }
        // The following test is just a heuristic, it is not required, and it is fine that is misses p == 0.
        if (mult != coeff_field_.additive_identity()) {  // For all columns in the boundary,
          for (auto cell_ref : col->col_) {  // insert every cell in map_a_ds with multiplicity
            Arith_element w_y = coeff_field_.times(cell_ref.coefficient_, mult);  // coefficient * multiplicity
  
            if (w_y != coeff_field_.additive_identity()) {  // if != 0
              result_insert_a_ds = map_a_ds.insert(std::pair<Simplex_key, Arith_element>(cell_ref.key_, w_y));
              if (!(result_insert_a_ds.second)) {  // if cell_ref.key_ already a Key in map_a_ds
                result_insert_a_ds.first->second = coeff_field_.plus_equal(result_insert_a_ds.first->second, w_y);
                if (result_insert_a_ds.first->second == coeff_field_.additive_identity()) {
                  map_a_ds.erase(result_insert_a_ds.first);
                }
              }
            }
          }
        }
      }
    }
  
    /*
     * Update the cohomology groups under the insertion of a simplex.
     */
    void update_cohomology_groups(Simplex_handle sigma, int dim_sigma) {
  // Compute the annotation of the boundary of sigma:
      std::map<Simplex_key, Arith_element> map_a_ds;
      annotation_of_the_boundary(map_a_ds, sigma, dim_sigma);
  // Update the cohomology groups:
      if (map_a_ds.empty()) {  // sigma is a creator in all fields represented in coeff_field_
        if (dim_sigma < dim_max_) {
          create_cocycle(sigma, coeff_field_.multiplicative_identity(),
                         coeff_field_.characteristic());
        }
      } else {        // sigma is a destructor in at least a field in coeff_field_
        // Convert map_a_ds to a vector
        A_ds_type a_ds;  // admits reverse iterators
        for (auto map_a_ds_ref : map_a_ds) {
          a_ds.push_back(
              std::pair<Simplex_key, Arith_element>(map_a_ds_ref.first,
                                                    map_a_ds_ref.second));
        }
  
        Arith_element inv_x, charac;
        Arith_element prod = coeff_field_.characteristic();  // Product of characteristic of the fields
        for (auto a_ds_rit = a_ds.rbegin();
            (a_ds_rit != a_ds.rend())
                && (prod != coeff_field_.multiplicative_identity()); ++a_ds_rit) {
          std::tie(inv_x, charac) = coeff_field_.inverse(a_ds_rit->second, prod);
  
          if (inv_x != coeff_field_.additive_identity()) {
            destroy_cocycle(sigma, a_ds, a_ds_rit->first, inv_x, charac);
            prod /= charac;
          }
        }
        if (prod != coeff_field_.multiplicative_identity()
            && dim_sigma < dim_max_) {
          create_cocycle(sigma, coeff_field_.multiplicative_identity(prod), prod);
        }
      }
    }
  
    /*  \brief Create a new cocycle class.
     *
     * The class is created by the insertion of the simplex sigma.
     * The methods adds a cocycle, representing the new cocycle class,
     * to the matrix representing the cohomology groups.
     * The new cocycle has value 0 on every simplex except on sigma
     * where it worths 1.*/
    void create_cocycle(Simplex_handle sigma, Arith_element x,
                        Arith_element charac) {
      Simplex_key key = cpx_->key(sigma);
      // Create a column containing only one cell,
      Column * new_col = column_pool_.construct(key);
      Cell * new_cell = cell_pool_.construct(key, x, new_col);
      new_col->col_.push_back(*new_cell);
      // and insert it in the matrix, in constant time thanks to the hint cam_.end().
      // Indeed *new_col has the biggest lexicographic value because key is the
      // biggest key used so far.
      cam_.insert(cam_.end(), *new_col);
      // Update the disjoint sets data structure.
      Hcell * new_hcell = new Hcell;
      new_hcell->push_back(*new_cell);
      transverse_idx_[key] = cocycle(charac, new_hcell);  // insert the new row
      ds_repr_[key] = new_col;
    }
  
    /*  \brief Destroy a cocycle class.
     *
     * The cocycle class is destroyed by the insertion of sigma.
     * The methods proceeds to a reduction of the matrix representing
     * the cohomology groups using Gauss pivoting. The reduction zeros-out
     * the row containing the cell with highest key in
     * a_ds, the annotation of the boundary of simplex sigma. This key
     * is "death_key".*/
    void destroy_cocycle(Simplex_handle sigma, A_ds_type const& a_ds,
                         Simplex_key death_key, Arith_element inv_x,
                         Arith_element charac) {
      // Create a finite persistent interval for which the interval exists
      if (interval_length_policy(cpx_->simplex(death_key), sigma)) {
        persistent_pairs_.emplace_back(cpx_->simplex(death_key)  // creator
            , sigma                                              // destructor
            , charac);                                           // fields
      }
  
      auto death_key_row = transverse_idx_.find(death_key);  // Find the beginning of the row.
      std::pair<typename Cam::iterator, bool> result_insert_cam;
  
      auto row_cell_it = death_key_row->second.row_->begin();
  
      while (row_cell_it != death_key_row->second.row_->end()) {  // Traverse all cells in
        // the row at index death_key.
        Arith_element w = coeff_field_.times_minus(inv_x, row_cell_it->coefficient_);
  
        if (w != coeff_field_.additive_identity()) {
          Column * curr_col = row_cell_it->self_col_;
          ++row_cell_it;
          // Disconnect the column from the rows in the CAM.
          for (auto& col_cell : curr_col->col_) {
            col_cell.base_hook_cam_h::unlink();
          }
  
          // Remove the column from the CAM before modifying its value
          cam_.erase(cam_.iterator_to(*curr_col));
          // Proceed to the reduction of the column
          plus_equal_column(*curr_col, a_ds, w);
  
          if (curr_col->col_.empty()) {  // If the column is null
            ds_repr_[curr_col->class_key_] = NULL;
            column_pool_.destroy(curr_col);  // delete curr_col;
          } else {
            // Find whether the column obtained is already in the CAM
            result_insert_cam = cam_.insert(*curr_col);
            if (result_insert_cam.second) {  // If it was not in the CAM before: insertion has succeeded
              for (auto& col_cell : curr_col->col_) {
                // re-establish the row links
                transverse_idx_[col_cell.key_].row_->push_front(col_cell);
              }
            } else {  // There is already an identical column in the CAM:
              // merge two disjoint sets.
              dsets_.link(curr_col->class_key_,
                          result_insert_cam.first->class_key_);
  
              Simplex_key key_tmp = dsets_.find_set(curr_col->class_key_);
              ds_repr_[key_tmp] = &(*(result_insert_cam.first));
              result_insert_cam.first->class_key_ = key_tmp;
              // intrusive containers don't own their elements, we have to release them manually
              curr_col->col_.clear_and_dispose([&](Cell*p){cell_pool_.destroy(p);});
              column_pool_.destroy(curr_col);  // delete curr_col;
            }
          }
        } else {
          ++row_cell_it;
        }  // If w == 0, pass.
      }
  
      // Because it is a killer simplex, set the data of sigma to null_key().
      if (charac == coeff_field_.characteristic()) {
        cpx_->assign_key(sigma, cpx_->null_key());
      }
      if (death_key_row->second.characteristics_ == charac) {
        delete death_key_row->second.row_;
        transverse_idx_.erase(death_key_row);
      } else {
        death_key_row->second.characteristics_ /= charac;
      }
    }
  
    /*
     * Assign:    target <- target + w * other.
     */
    void plus_equal_column(Column & target, A_ds_type const& other  // value_type is pair<Simplex_key,Arith_element>
                           , Arith_element w) {
      auto target_it = target.col_.begin();
      auto other_it = other.begin();
      while (target_it != target.col_.end() && other_it != other.end()) {
        if (target_it->key_ < other_it->first) {
          ++target_it;
        } else {
          if (target_it->key_ > other_it->first) {
            Cell * cell_tmp = cell_pool_.construct(Cell(other_it->first   // key
                , coeff_field_.additive_identity(), &target));
  
            cell_tmp->coefficient_ = coeff_field_.plus_times_equal(cell_tmp->coefficient_, other_it->second, w);
  
            target.col_.insert(target_it, *cell_tmp);
  
            ++other_it;
          } else {  // it1->key == it2->key
            // target_it->coefficient_ <- target_it->coefficient_ + other_it->second * w
            target_it->coefficient_ = coeff_field_.plus_times_equal(target_it->coefficient_, other_it->second, w);
            if (target_it->coefficient_ == coeff_field_.additive_identity()) {
              auto tmp_it = target_it;
              ++target_it;
              ++other_it;   // iterators remain valid
              Cell * tmp_cell_ptr = &(*tmp_it);
              target.col_.erase(tmp_it);  // removed from column
  
              cell_pool_.destroy(tmp_cell_ptr);  // delete from memory
            } else {
              ++target_it;
              ++other_it;
            }
          }
        }
      }
      while (other_it != other.end()) {
        Cell * cell_tmp = cell_pool_.construct(Cell(other_it->first, coeff_field_.additive_identity(), &target));
        cell_tmp->coefficient_ = coeff_field_.plus_times_equal(cell_tmp->coefficient_, other_it->second, w);
        target.col_.insert(target.col_.end(), *cell_tmp);
  
        ++other_it;
      }
    }
  
    /*
     * Compare two intervals by length.
     */
    struct cmp_intervals_by_length {
      explicit cmp_intervals_by_length(FilteredComplex * sc)
          : sc_(sc) {
      }
      bool operator()(const Persistent_interval & p1, const Persistent_interval & p2) {
        return (sc_->filtration(get < 1 > (p1)) - sc_->filtration(get < 0 > (p1))
            > sc_->filtration(get < 1 > (p2)) - sc_->filtration(get < 0 > (p2)));
      }
      FilteredComplex * sc_;
    };
  
   public:
    /** \brief Output the persistence diagram in ostream.
     *
     * The file format is the following:
     *    p1*...*pr   dim b d
     *
     * where "dim" is the dimension of the homological feature,
     * b and d are respectively the birth and death of the feature and
     * p1*...*pr is the product of prime numbers pi such that the homology
     * feature exists in homology with Z/piZ coefficients.
     */
    void output_diagram(std::ostream& ostream = std::cout) {
      cmp_intervals_by_length cmp(cpx_);
      std::sort(std::begin(persistent_pairs_), std::end(persistent_pairs_), cmp);
      bool has_infinity = std::numeric_limits<Filtration_value>::has_infinity;
      for (auto pair : persistent_pairs_) {
        // Special case on windows, inf is "1.#INF" (cf. unitary tests and R package TDA)
        if (has_infinity && cpx_->filtration(get<1>(pair)) == std::numeric_limits<Filtration_value>::infinity()) {
          ostream << get<2>(pair) << "  " << cpx_->dimension(get<0>(pair)) << " "
            << cpx_->filtration(get<0>(pair)) << " inf " << std::endl;
        } else {
          ostream << get<2>(pair) << "  " << cpx_->dimension(get<0>(pair)) << " "
            << cpx_->filtration(get<0>(pair)) << " "
            << cpx_->filtration(get<1>(pair)) << " " << std::endl;
        }
      }
    }
  
    void write_output_diagram(std::string diagram_name) {
      std::ofstream diagram_out(diagram_name.c_str());
      cmp_intervals_by_length cmp(cpx_);
      std::sort(std::begin(persistent_pairs_), std::end(persistent_pairs_), cmp);
      bool has_infinity = std::numeric_limits<Filtration_value>::has_infinity;
      for (auto pair : persistent_pairs_) {
        // Special case on windows, inf is "1.#INF"
        if (has_infinity && cpx_->filtration(get<1>(pair)) == std::numeric_limits<Filtration_value>::infinity()) {
          diagram_out << cpx_->dimension(get<0>(pair)) << " "
                << cpx_->filtration(get<0>(pair)) << " inf" << std::endl;
        } else {
          diagram_out << cpx_->dimension(get<0>(pair)) << " "
                << cpx_->filtration(get<0>(pair)) << " "
                << cpx_->filtration(get<1>(pair)) << std::endl;
        }
      }
    }
  
    /** @brief Returns Betti numbers.
     * @return A vector of Betti numbers.
     */
    std::vector<int> betti_numbers() const {
      // Don't allocate a vector of negative size for an empty complex
      int siz = std::max(dim_max_, 0);
      // Init Betti numbers vector with zeros until Simplicial complex dimension
      std::vector<int> betti_numbers(siz);
  
      for (auto pair : persistent_pairs_) {
        // Count never ended persistence intervals
        if (cpx_->null_simplex() == get<1>(pair)) {
          // Increment corresponding betti number
          betti_numbers[cpx_->dimension(get<0>(pair))] += 1;
        }
      }
      return betti_numbers;
    }
  
    /** @brief Returns the Betti number of the dimension passed by parameter.
     * @param[in] dimension The Betti number dimension to get.
     * @return Betti number of the given dimension
     *
     */
    int betti_number(int dimension) const {
      int betti_number = 0;
  
      for (auto pair : persistent_pairs_) {
        // Count never ended persistence intervals
        if (cpx_->null_simplex() == get<1>(pair)) {
          if (cpx_->dimension(get<0>(pair)) == dimension) {
            // Increment betti number found
            ++betti_number;
          }
        }
      }
      return betti_number;
    }
  
    /** @brief Returns the persistent Betti numbers.
     * @param[in] from The persistence birth limit to be added in the number \f$(persistent birth \leq from)\f$.
     * @param[in] to The persistence death limit to be added in the number  \f$(persistent death > to)\f$.
     * @return A vector of persistent Betti numbers.
     */
    std::vector<int> persistent_betti_numbers(Filtration_value from, Filtration_value to) const {
      // Don't allocate a vector of negative size for an empty complex
      int siz = std::max(dim_max_, 0);
      // Init Betti numbers vector with zeros until Simplicial complex dimension
      std::vector<int> betti_numbers(siz);
      for (auto pair : persistent_pairs_) {
        // Count persistence intervals that covers the given interval
        // null_simplex test : if the function is called with to=+infinity, we still get something useful. And it will
        // still work if we change the complex filtration function to reject null simplices.
        if (cpx_->filtration(get<0>(pair)) <= from &&
            (get<1>(pair) == cpx_->null_simplex() || cpx_->filtration(get<1>(pair)) > to)) {
          // Increment corresponding betti number
          betti_numbers[cpx_->dimension(get<0>(pair))] += 1;
        }
      }
      return betti_numbers;
    }
  
    /** @brief Returns the persistent Betti number of the dimension passed by parameter.
     * @param[in] dimension The Betti number dimension to get.
     * @param[in] from The persistence birth limit to be added in the number \f$(persistent birth \leq from)\f$.
     * @param[in] to The persistence death limit to be added in the number  \f$(persistent death > to)\f$.
     * @return Persistent Betti number of the given dimension
     */
    int persistent_betti_number(int dimension, Filtration_value from, Filtration_value to) const {
      int betti_number = 0;
  
      for (auto pair : persistent_pairs_) {
        // Count persistence intervals that covers the given interval
        // null_simplex test : if the function is called with to=+infinity, we still get something useful. And it will
        // still work if we change the complex filtration function to reject null simplices.
        if (cpx_->filtration(get<0>(pair)) <= from &&
            (get<1>(pair) == cpx_->null_simplex() || cpx_->filtration(get<1>(pair)) > to)) {
          if (cpx_->dimension(get<0>(pair)) == dimension) {
            // Increment betti number found
            ++betti_number;
          }
        }
      }
      return betti_number;
    }
  
    /** @brief Returns a list of persistence birth and death FilteredComplex::Simplex_handle pairs.
     * @return A list of Persistent_cohomology::Persistent_interval
     */
    const std::vector<Persistent_interval>& get_persistent_pairs() const {
      return persistent_pairs_;
    }
  
    /** @brief Returns persistence intervals for a given dimension.
     * @param[in] dimension Dimension to get the birth and death pairs from.
     * @return A vector of persistence intervals (birth and death) on a fixed dimension.
     */
    std::vector< std::pair< Filtration_value , Filtration_value > >
    intervals_in_dimension(int dimension) {
      std::vector< std::pair< Filtration_value , Filtration_value > > result;
      // auto && pair, to avoid unnecessary copying
      for (auto && pair : persistent_pairs_) {
        if (cpx_->dimension(get<0>(pair)) == dimension) {
          result.emplace_back(cpx_->filtration(get<0>(pair)), cpx_->filtration(get<1>(pair)));
        }
      }
      return result;
    }
  
   private:
    /*
     * Structure representing a cocycle.
     */
    struct cocycle {
      cocycle()
          : row_(nullptr),
            characteristics_() {
      }
      cocycle(Arith_element characteristics, Hcell * row)
          : row_(row),
            characteristics_(characteristics) {
      }
  
      Hcell * row_;                    // points to the corresponding row in the CAM
      Arith_element characteristics_;  // product of field characteristics for which the cocycle exist
    };
  
   public:
    FilteredComplex * cpx_;
    int dim_max_;
    CoefficientField coeff_field_;
    size_t num_simplices_;
  
    /*  Disjoint sets data structure to link the model of FilteredComplex
     * with the compressed annotation matrix.
     * ds_rank_ is a property map Simplex_key -> int, ds_parent_ is a property map
     * Simplex_key -> simplex_key_t */
    std::vector<int> ds_rank_;
    std::vector<Simplex_key> ds_parent_;
    std::vector<Column *> ds_repr_;
    boost::disjoint_sets<int *, Simplex_key *> dsets_;
    /* The compressed annotation matrix fields.*/
    Cam cam_;
    /*  Dictionary establishing the correspondance between the Simplex_key of
     * the root vertex in the union-find ds and the Simplex_key of the vertex which
     * created the connected component as a 0-dimension homology feature.*/
    std::map<Simplex_key, Simplex_key> zero_cocycles_;
    /*  Key -> row. */
    std::map<Simplex_key, cocycle> transverse_idx_;
    /* Persistent intervals. */
    std::vector<Persistent_interval> persistent_pairs_;
    length_interval interval_length_policy;
  
    Simple_object_pool<Column> column_pool_;
    Simple_object_pool<Cell> cell_pool_;
  };
  
  }  // namespace persistent_cohomology
  
  }  // namespace Gudhi
  
  #endif  // PERSISTENT_COHOMOLOGY_H_
  