/*    This file is part of the Gudhi Library. The Gudhi library
 *    (Geometric Understanding in Higher Dimensions) is a generic C++
 *    library for computational topology.
 *
 *    Author(s):       Clément Maria
 *
 *    Copyright (C) 2014  INRIA Sophia Antipolis-Méditerranée (France)
 *
 *    This program is free software: you can redistribute it and/or modify
 *    it under the terms of the GNU General Public License as published by
 *    the Free Software Foundation, either version 3 of the License, or
 *    (at your option) any later version.
 *
 *    This program is distributed in the hope that it will be useful,
 *    but WITHOUT ANY WARRANTY; without even the implied warranty of
 *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *    GNU General Public License for more details.
 *
 *    You should have received a copy of the GNU General Public License
 *    along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */

#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>

namespace Gudhi {

namespace persistent_cohomology {

/** \defgroup persistent_cohomology Persistent Cohomology
  *

  Computation of persistent cohomology using the algorithm of 
   \cite DBLP:journals/dcg/SilvaMV11 and \cite DBLP:journals/corr/abs-1208-5018 
   and the Compressed Annotation Matrix 
   implementation of \cite DBLP:conf/esa/BoissonnatDM13 
       
  The theory of homology consists in attaching to a topological space a sequence of 
  (homology) groups, 
  capturing global topological features 
  like connected components, holes, cavities, etc. Persistent homology studies the evolution 
  -- birth, life and death -- of 
  these features when the topological space is changing. Consequently, the theory is essentially 
  composed of three elements: 
  topological spaces, their homology groups and an evolution scheme.

 The theory of homology consists in attaching to a topological space a sequence of
 (homology) groups,
 capturing global topological features
 like connected components, holes, cavities, etc. Persistent homology studies the evolution
 -- birth, life and death -- of
 these features when the topological space is changing. Consequently, the theory is essentially
 composed of three elements:
 topological spaces, their homology groups and an evolution scheme.

 <DT>Topological Spaces:</DT>
 Topological spaces are represented by simplicial complexes.
 Let \f$V = \{1, \cdots ,|V|\}\f$ be a set of <EM>vertices</EM>.
 A <EM>simplex</EM> \f$\sigma\f$ is a subset of vertices
 \f$\sigma \subseteq V\f$. A <EM>simplicial complex</EM> \f$\mathbf{K}\f$
 on \f$V\f$ is a collection of simplices \f$\{\sigma\}\f$,
 \f$\sigma \subseteq V\f$, such that \f$\tau \subseteq \sigma \in \mathbf{K}
 \Rightarrow \tau \in \mathbf{K}\f$. The dimension \f$n=|\sigma|-1\f$ of \f$\sigma\f$
 is its number of elements minus 1. A <EM>filtration</EM> of a simplicial complex is
 a function \f$f:\mathbf{K} \rightarrow \mathbb{R}\f$ satisfying \f$f(\tau)\leq
 f(\sigma)\f$ whenever \f$\tau \subseteq \sigma\f$.

 We define the concept FilteredComplex which enumerates the requirements for a class
 to represent a filtered complex from which persistent homology may be computed.
 We use the vocabulary of simplicial complexes, but the concept
 is valid for any type of cell complex. The main requirements
 are the definition of:
 \li type <CODE>Indexing_tag</CODE>, which is a model of the concept
 <CODE>IndexingTag</CODE>,
 describing the nature of the indexing scheme,
 \li type Simplex_handle to manipulate simplices,
 \li method <CODE>int dimension(Simplex_handle)</CODE> returning
 the dimension of a simplex,
 \li type and method <CODE>Boundary_simplex_range
 boundary_simplex_range(Simplex_handle)</CODE> that returns
 a range giving access to the codimension 1 subsimplices of the
 input simplex, as-well-as the coefficients \f$(-1)^i\f$ in the
 definition of the operator \f$\partial\f$. The iterators have
 value type <CODE>Simplex_handle</CODE>,
 \li type and method
 <CODE>Filtration_simplex_range filtration_simplex_range ()</CODE>
 that returns a range giving
 access to all the simplices of the complex read in the order
 assigned by the indexing scheme,
 \li type and method
 <CODE>Filtration_value filtration (Simplex_handle)</CODE> that returns the value of
 the filtration on the simplex represented by the handle.

 <DT>Homology:</DT>
 For a ring \f$\mathcal{R}\f$, the group of <EM>n-chains</EM>,
 denoted \f$\mathbf{C}_n(\mathbf{K},\mathcal{R})\f$, of \f$\mathbf{K}\f$ is the
 group of formal sums of
 n-simplices with \f$\mathcal{R}\f$ coefficients. The <EM>boundary operator</EM> is a
 linear operator
 \f$\partial_n: \mathbf{C}_n(\mathbf{K},\mathcal{R}) \rightarrow \mathbf{C}_{n-1}(\mathbf{K},\mathcal{R})\f$
 such that \f$\partial_n \sigma = \partial_n [v_0, \cdots , v_n] =
 \sum_{i=0}^n (-1)^{i}[v_0,\cdots ,\widehat{v_i}, \cdots,v_n]\f$,
 where \f$\widehat{v_i}\f$ means \f$v_i\f$ is omitted from the list. The chain
 groups form a sequence:

 \f[\cdots \ \ \mathbf{C}_n(\mathbf{K},\mathcal{R}) \xrightarrow{\ \partial_n\ } \mathbf{C}_{n-1}(\mathbf{K},\mathcal{R})
 \xrightarrow{\partial_{n-1}} \cdots \xrightarrow{\ \partial_2 \ }
 \mathbf{C}_1(\mathbf{K},\mathcal{R}) \xrightarrow{\ \partial_1 \ }  \mathbf{C}_0(\mathbf{K},\mathcal{R}) \f]

 of finitely many groups \f$\mathbf{C}_n(\mathbf{K},\mathcal{R})\f$ and homomorphisms
 \f$\partial_n\f$, indexed by the dimension \f$n \geq 0\f$.
 The boundary operators satisfy the property \f$\partial_n \circ \partial_{n+1}=0\f$
 for every \f$n > 0\f$
 and we define the homology groups:

 \f[\mathbf{H}_n(\mathbf{K},\mathcal{R}) = \ker \partial_n / \mathrm{im} \  \partial_{n+1}\f]

 We refer to \cite Munkres-elementsalgtop1984 for an introduction to homology
 theory and to \cite DBLP:books/daglib/0025666 for an introduction to persistent homology.

 <DT>Indexing Scheme:</DT>
 "Changing" a simplicial complex consists in applying a simplicial map.
 An <EM>indexing scheme</EM> is a directed graph together with a traversal
 order, such that two
 consecutive nodes in the graph are connected by an arrow (either forward or backward).
 The nodes represent simplicial complexes and the directed edges simplicial maps.

 From the computational point of view, there are two types of indexing schemes of
 interest
 in persistent homology: <EM>linear</EM> ones
 \f$\bullet \longrightarrow \bullet \longrightarrow \cdots \longrightarrow \bullet
 \longrightarrow \bullet\f$
 in persistent homology \cite DBLP:journals/dcg/ZomorodianC05 ,
 and <EM>zigzag</EM> ones
 \f$\bullet \longrightarrow \bullet \longleftarrow \cdots
 \longrightarrow \bullet
 \longleftarrow \bullet \f$ in zigzag persistent
 homology \cite DBLP:journals/focm/CarlssonS10.
 These indexing schemes have a natural left-to-right traversal order, and we
 describe them with ranges and iterators.
 In the current release of the Gudhi library, only the linear case is implemented.

 In the following, we consider the case where the indexing scheme is induced
 by a filtration.
 Ordering the simplices
 by increasing filtration values (breaking ties so as a simplex appears after
 its subsimplices of same filtration value) provides an indexing scheme.

\section Examples
    We provide several example files: run these examples with -h for details on their use, and read the README file.

\li <CODE>rips_persistence.cpp</CODE> computes the Rips complex of a point cloud and its persistence diagram.

\li <CODE>rips_multifield_persistence.cpp</CODE> computes the Rips complex of a point cloud and its persistence diagram 
with a family of field coefficients.

\li <CODE>performance_rips_persistence.cpp</CODE> provides timings for the construction of the Rips complex on a set of 
points sampling a Klein bottle in \f$\mathbb{R}^5\f$ with a simplex tree, its conversion to a 
Hasse diagram and the computation of persistent homology and multi-field persistent homology for the 
different representations.



 \author    Clément Maria
 \version   1.0
 \date      2014
 \copyright GNU General Public License v3.
 @{  
 */

/** \brief Computes the persistent cohomology of a filtered complex.
 *
 * 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
 *
 */
// Memory allocation policy: classic, use a mempool, etc.*/
template<class FilteredComplex, class CoefficientField>  // to do mem allocation policy: classic, mempool, etc.
class Persistent_cohomology {
 public:
  typedef FilteredComplex Complex_ds;
  // Data attached to each simplex to interface with a Property Map.
  typedef typename Complex_ds::Simplex_key Simplex_key;
  typedef typename Complex_ds::Simplex_handle Simplex_handle;
  typedef typename Complex_ds::Filtration_value Filtration_value;
  typedef typename CoefficientField::Element Arith_element;
// 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;
// Persistent interval type. The Arith_element field is used for the multi-field framework.
  typedef std::tuple<Simplex_handle, Simplex_handle, Arith_element> Persistent_interval;

  /** \brief Initializes the Persistent_cohomology class.
   *
   * @param[in] cpx Complex for which the persistent homology is computed.
   cpx is a model of FilteredComplex
   */
  explicit Persistent_cohomology(Complex_ds& cpx)
      : 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_() {
    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));
    }
  }

  /** \brief Initializes the Persistent_cohomology class.
   *
   * @param[in] cpx Complex for which the persistent homology is compiuted.
   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.
   */
  Persistent_cohomology(Complex_ds& cpx, bool persistence_dim_max)
      : Persistent_cohomology(cpx) {
    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(Complex_ds * 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;
    }

    Complex_ds * 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 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;
    // Danger: not thread-safe!
    static 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(Complex_ds * 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)));
    }
    Complex_ds * 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);
    for (auto pair : persistent_pairs_) {
    diagram_out << cpx_->dimension(get<0>(pair)) << " "
          << cpx_->filtration(get<0>(pair)) << " "
          << cpx_->filtration(get<1>(pair)) << std::endl;
    }
  }

 private:
  /*
   * Structure representing a cocycle.
   */
  struct cocycle {
    cocycle() {
    }
    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:
  Complex_ds * 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_;
};

/** @} */  // end defgroup persistent_cohomology

}  // namespace persistent_cohomology

}  // namespace Gudhi

#endif  // PERSISTENT_COHOMOLOGY_H_