Persistent_cohomology.h
1 /* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT.
2  * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details.
3  * Author(s): ClĂ©ment Maria
4  *
5  * Copyright (C) 2014 Inria
6  *
7  * Modification(s):
8  * - YYYY/MM Author: Description of the modification
9  */
10 
11 #ifndef PERSISTENT_COHOMOLOGY_H_
12 #define PERSISTENT_COHOMOLOGY_H_
13 
14 #include <gudhi/Persistent_cohomology/Persistent_cohomology_column.h>
15 #include <gudhi/Persistent_cohomology/Field_Zp.h>
16 #include <gudhi/Simple_object_pool.h>
17 
18 #include <boost/intrusive/set.hpp>
19 #include <boost/pending/disjoint_sets.hpp>
20 #include <boost/intrusive/list.hpp>
21 
22 #include <map>
23 #include <utility>
24 #include <list>
25 #include <vector>
26 #include <set>
27 #include <fstream> // std::ofstream
28 #include <limits> // for numeric_limits<>
29 #include <tuple>
30 #include <algorithm>
31 #include <string>
32 #include <stdexcept> // for std::out_of_range
33 
34 namespace Gudhi {
35 
36 namespace persistent_cohomology {
37 
50 // TODO(CM): Memory allocation policy: classic, use a mempool, etc.
51 template<class FilteredComplex, class CoefficientField>
53  public:
54  // Data attached to each simplex to interface with a Property Map.
55 
66  typedef std::tuple<Simplex_handle, Simplex_handle, Arith_element> Persistent_interval;
67 
68  private:
69  // Compressed Annotation Matrix types:
70  // Column type
71  typedef Persistent_cohomology_column<Simplex_key, Arith_element> Column; // contains 1 set_hook
72  // Cell type
73  typedef typename Column::Cell Cell; // contains 2 list_hooks
74  // Remark: constant_time_size must be false because base_hook_cam_h has auto_unlink link_mode
75  typedef boost::intrusive::list<Cell,
76  boost::intrusive::constant_time_size<false>,
77  boost::intrusive::base_hook<base_hook_cam_h> > Hcell;
78 
79  typedef boost::intrusive::set<Column,
80  boost::intrusive::constant_time_size<false> > Cam;
81  // Sparse column type for the annotation of the boundary of an element.
82  typedef std::vector<std::pair<Simplex_key, Arith_element> > A_ds_type;
83 
84  public:
95  explicit Persistent_cohomology(FilteredComplex& cpx, bool persistence_dim_max = false)
96  : cpx_(&cpx),
97  dim_max_(cpx.dimension()), // upper bound on the dimension of the simplices
98  coeff_field_(), // initialize the field coefficient structure.
99  num_simplices_(cpx_->num_simplices()), // num_simplices save to avoid to call thrice the function
100  ds_rank_(num_simplices_), // union-find
101  ds_parent_(num_simplices_), // union-find
102  ds_repr_(num_simplices_, NULL), // union-find -> annotation vectors
103  dsets_(&ds_rank_[0], &ds_parent_[0]), // union-find
104  cam_(), // collection of annotation vectors
105  zero_cocycles_(), // union-find -> Simplex_key of creator for 0-homology
106  transverse_idx_(), // key -> row
107  persistent_pairs_(),
108  interval_length_policy(&cpx, 0),
109  column_pool_(), // memory pools for the CAM
110  cell_pool_() {
111  if (cpx_->num_simplices() > std::numeric_limits<Simplex_key>::max()) {
112  // num_simplices must be strictly lower than the limit, because a value is reserved for null_key.
113  throw std::out_of_range("The number of simplices is more than Simplex_key type numeric limit.");
114  }
115  Simplex_key idx_fil = 0;
116  for (auto sh : cpx_->filtration_simplex_range()) {
117  cpx_->assign_key(sh, idx_fil);
118  ++idx_fil;
119  dsets_.make_set(cpx_->key(sh));
120  }
121  if (persistence_dim_max) {
122  ++dim_max_;
123  }
124  }
125 
127  // Clean the transversal lists
128  for (auto & transverse_ref : transverse_idx_) {
129  // Destruct all the cells
130  transverse_ref.second.row_->clear_and_dispose([&](Cell*p){p->~Cell();});
131  delete transverse_ref.second.row_;
132  }
133  }
134 
135  private:
136  struct length_interval {
137  length_interval(FilteredComplex * cpx, Filtration_value min_length)
138  : cpx_(cpx),
139  min_length_(min_length) {
140  }
141 
142  bool operator()(Simplex_handle sh1, Simplex_handle sh2) {
143  return cpx_->filtration(sh2) - cpx_->filtration(sh1) > min_length_;
144  }
145 
146  void set_length(Filtration_value new_length) {
147  min_length_ = new_length;
148  }
149 
150  FilteredComplex * cpx_;
151  Filtration_value min_length_;
152  };
153 
154  public:
156  void init_coefficients(int charac) {
157  coeff_field_.init(charac);
158  }
160  void init_coefficients(int charac_min, int charac_max) {
161  coeff_field_.init(charac_min, charac_max);
162  }
163 
172  void compute_persistent_cohomology(Filtration_value min_interval_length = 0) {
173  interval_length_policy.set_length(min_interval_length);
174  // Compute all finite intervals
175  for (auto sh : cpx_->filtration_simplex_range()) {
176  int dim_simplex = cpx_->dimension(sh);
177  switch (dim_simplex) {
178  case 0:
179  break;
180  case 1:
181  update_cohomology_groups_edge(sh);
182  break;
183  default:
184  update_cohomology_groups(sh, dim_simplex);
185  break;
186  }
187  }
188  // Compute infinite intervals of dimension 0
189  Simplex_key key;
190  for (auto v_sh : cpx_->skeleton_simplex_range(0)) { // for all 0-dimensional simplices
191  key = cpx_->key(v_sh);
192 
193  if (ds_parent_[key] == key // root of its tree
194  && zero_cocycles_.find(key) == zero_cocycles_.end()) {
195  persistent_pairs_.emplace_back(
196  cpx_->simplex(key), cpx_->null_simplex(), coeff_field_.characteristic());
197  }
198  }
199  for (auto zero_idx : zero_cocycles_) {
200  persistent_pairs_.emplace_back(
201  cpx_->simplex(zero_idx.second), cpx_->null_simplex(), coeff_field_.characteristic());
202  }
203  // Compute infinite interval of dimension > 0
204  for (auto cocycle : transverse_idx_) {
205  persistent_pairs_.emplace_back(
206  cpx_->simplex(cocycle.first), cpx_->null_simplex(), cocycle.second.characteristics_);
207  }
208  }
209 
210  private:
215  void update_cohomology_groups_edge(Simplex_handle sigma) {
216  Simplex_handle u, v;
217  boost::tie(u, v) = cpx_->endpoints(sigma);
218 
219  Simplex_key ku = dsets_.find_set(cpx_->key(u));
220  Simplex_key kv = dsets_.find_set(cpx_->key(v));
221 
222  if (ku != kv) { // Destroy a connected component
223  dsets_.link(ku, kv);
224  // Keys of the simplices which created the connected components containing
225  // respectively u and v.
226  Simplex_key idx_coc_u, idx_coc_v;
227  auto map_it_u = zero_cocycles_.find(ku);
228  // If the index of the cocycle representing the class is already ku.
229  if (map_it_u == zero_cocycles_.end()) {
230  idx_coc_u = ku;
231  } else {
232  idx_coc_u = map_it_u->second;
233  }
234 
235  auto map_it_v = zero_cocycles_.find(kv);
236  // If the index of the cocycle representing the class is already kv.
237  if (map_it_v == zero_cocycles_.end()) {
238  idx_coc_v = kv;
239  } else {
240  idx_coc_v = map_it_v->second;
241  }
242 
243  if (cpx_->filtration(cpx_->simplex(idx_coc_u))
244  < cpx_->filtration(cpx_->simplex(idx_coc_v))) { // Kill cocycle [idx_coc_v], which is younger.
245  if (interval_length_policy(cpx_->simplex(idx_coc_v), sigma)) {
246  persistent_pairs_.emplace_back(
247  cpx_->simplex(idx_coc_v), sigma, coeff_field_.characteristic());
248  }
249  // Maintain the index of the 0-cocycle alive.
250  if (kv != idx_coc_v) {
251  zero_cocycles_.erase(map_it_v);
252  }
253  if (kv == dsets_.find_set(kv)) {
254  if (ku != idx_coc_u) {
255  zero_cocycles_.erase(map_it_u);
256  }
257  zero_cocycles_[kv] = idx_coc_u;
258  }
259  } else { // Kill cocycle [idx_coc_u], which is younger.
260  if (interval_length_policy(cpx_->simplex(idx_coc_u), sigma)) {
261  persistent_pairs_.emplace_back(
262  cpx_->simplex(idx_coc_u), sigma, coeff_field_.characteristic());
263  }
264  // Maintain the index of the 0-cocycle alive.
265  if (ku != idx_coc_u) {
266  zero_cocycles_.erase(map_it_u);
267  }
268  if (ku == dsets_.find_set(ku)) {
269  if (kv != idx_coc_v) {
270  zero_cocycles_.erase(map_it_v);
271  }
272  zero_cocycles_[ku] = idx_coc_v;
273  }
274  }
275  cpx_->assign_key(sigma, cpx_->null_key());
276  } else if (dim_max_ > 1) { // If ku == kv, same connected component: create a 1-cocycle class.
277  create_cocycle(sigma, coeff_field_.multiplicative_identity(), coeff_field_.characteristic());
278  }
279  }
280 
281  /*
282  * Compute the annotation of the boundary of a simplex.
283  */
284  void annotation_of_the_boundary(
285  std::map<Simplex_key, Arith_element> & map_a_ds, Simplex_handle sigma,
286  int dim_sigma) {
287  // traverses the boundary of sigma, keeps track of the annotation vectors,
288  // with multiplicity. We used to sum the coefficients directly in
289  // annotations_in_boundary by using a map, we now do it later.
290  typedef std::pair<Column *, int> annotation_t;
291 #ifdef GUDHI_CAN_USE_CXX11_THREAD_LOCAL
292  thread_local
293 #endif // GUDHI_CAN_USE_CXX11_THREAD_LOCAL
294  std::vector<annotation_t> annotations_in_boundary;
295  annotations_in_boundary.clear();
296  int sign = 1 - 2 * (dim_sigma % 2); // \in {-1,1} provides the sign in the
297  // alternate sum in the boundary.
298  Simplex_key key;
299  Column * curr_col;
300 
301  for (auto sh : cpx_->boundary_simplex_range(sigma)) {
302  key = cpx_->key(sh);
303  if (key != cpx_->null_key()) { // A simplex with null_key is a killer, and have null annotation
304  // Find its annotation vector
305  curr_col = ds_repr_[dsets_.find_set(key)];
306  if (curr_col != NULL) { // and insert it in annotations_in_boundary with multyiplicative factor "sign".
307  annotations_in_boundary.emplace_back(curr_col, sign);
308  }
309  }
310  sign = -sign;
311  }
312  // Place identical annotations consecutively so we can easily sum their multiplicities.
313  std::sort(annotations_in_boundary.begin(), annotations_in_boundary.end(),
314  [](annotation_t const& a, annotation_t const& b) { return a.first < b.first; });
315 
316  // Sum the annotations with multiplicity, using a map<key,coeff>
317  // to represent a sparse vector.
318  std::pair<typename std::map<Simplex_key, Arith_element>::iterator, bool> result_insert_a_ds;
319 
320  for (auto ann_it = annotations_in_boundary.begin(); ann_it != annotations_in_boundary.end(); ) {
321  Column* col = ann_it->first;
322  int mult = ann_it->second;
323  while (++ann_it != annotations_in_boundary.end() && ann_it->first == col) {
324  mult += ann_it->second;
325  }
326  // The following test is just a heuristic, it is not required, and it is fine that is misses p == 0.
327  if (mult != coeff_field_.additive_identity()) { // For all columns in the boundary,
328  for (auto cell_ref : col->col_) { // insert every cell in map_a_ds with multiplicity
329  Arith_element w_y = coeff_field_.times(cell_ref.coefficient_, mult); // coefficient * multiplicity
330 
331  if (w_y != coeff_field_.additive_identity()) { // if != 0
332  result_insert_a_ds = map_a_ds.insert(std::pair<Simplex_key, Arith_element>(cell_ref.key_, w_y));
333  if (!(result_insert_a_ds.second)) { // if cell_ref.key_ already a Key in map_a_ds
334  result_insert_a_ds.first->second = coeff_field_.plus_equal(result_insert_a_ds.first->second, w_y);
335  if (result_insert_a_ds.first->second == coeff_field_.additive_identity()) {
336  map_a_ds.erase(result_insert_a_ds.first);
337  }
338  }
339  }
340  }
341  }
342  }
343  }
344 
345  /*
346  * Update the cohomology groups under the insertion of a simplex.
347  */
348  void update_cohomology_groups(Simplex_handle sigma, int dim_sigma) {
349 // Compute the annotation of the boundary of sigma:
350  std::map<Simplex_key, Arith_element> map_a_ds;
351  annotation_of_the_boundary(map_a_ds, sigma, dim_sigma);
352 // Update the cohomology groups:
353  if (map_a_ds.empty()) { // sigma is a creator in all fields represented in coeff_field_
354  if (dim_sigma < dim_max_) {
355  create_cocycle(sigma, coeff_field_.multiplicative_identity(),
356  coeff_field_.characteristic());
357  }
358  } else { // sigma is a destructor in at least a field in coeff_field_
359  // Convert map_a_ds to a vector
360  A_ds_type a_ds; // admits reverse iterators
361  for (auto map_a_ds_ref : map_a_ds) {
362  a_ds.push_back(
363  std::pair<Simplex_key, Arith_element>(map_a_ds_ref.first,
364  map_a_ds_ref.second));
365  }
366 
367  Arith_element inv_x, charac;
368  Arith_element prod = coeff_field_.characteristic(); // Product of characteristic of the fields
369  for (auto a_ds_rit = a_ds.rbegin();
370  (a_ds_rit != a_ds.rend())
371  && (prod != coeff_field_.multiplicative_identity()); ++a_ds_rit) {
372  std::tie(inv_x, charac) = coeff_field_.inverse(a_ds_rit->second, prod);
373 
374  if (inv_x != coeff_field_.additive_identity()) {
375  destroy_cocycle(sigma, a_ds, a_ds_rit->first, inv_x, charac);
376  prod /= charac;
377  }
378  }
379  if (prod != coeff_field_.multiplicative_identity()
380  && dim_sigma < dim_max_) {
381  create_cocycle(sigma, coeff_field_.multiplicative_identity(prod), prod);
382  }
383  }
384  }
385 
386  /* \brief Create a new cocycle class.
387  *
388  * The class is created by the insertion of the simplex sigma.
389  * The methods adds a cocycle, representing the new cocycle class,
390  * to the matrix representing the cohomology groups.
391  * The new cocycle has value 0 on every simplex except on sigma
392  * where it worths 1.*/
393  void create_cocycle(Simplex_handle sigma, Arith_element x,
394  Arith_element charac) {
395  Simplex_key key = cpx_->key(sigma);
396  // Create a column containing only one cell,
397  Column * new_col = column_pool_.construct(key);
398  Cell * new_cell = cell_pool_.construct(key, x, new_col);
399  new_col->col_.push_back(*new_cell);
400  // and insert it in the matrix, in constant time thanks to the hint cam_.end().
401  // Indeed *new_col has the biggest lexicographic value because key is the
402  // biggest key used so far.
403  cam_.insert(cam_.end(), *new_col);
404  // Update the disjoint sets data structure.
405  Hcell * new_hcell = new Hcell;
406  new_hcell->push_back(*new_cell);
407  transverse_idx_[key] = cocycle(charac, new_hcell); // insert the new row
408  ds_repr_[key] = new_col;
409  }
410 
411  /* \brief Destroy a cocycle class.
412  *
413  * The cocycle class is destroyed by the insertion of sigma.
414  * The methods proceeds to a reduction of the matrix representing
415  * the cohomology groups using Gauss pivoting. The reduction zeros-out
416  * the row containing the cell with highest key in
417  * a_ds, the annotation of the boundary of simplex sigma. This key
418  * is "death_key".*/
419  void destroy_cocycle(Simplex_handle sigma, A_ds_type const& a_ds,
420  Simplex_key death_key, Arith_element inv_x,
421  Arith_element charac) {
422  // Create a finite persistent interval for which the interval exists
423  if (interval_length_policy(cpx_->simplex(death_key), sigma)) {
424  persistent_pairs_.emplace_back(cpx_->simplex(death_key) // creator
425  , sigma // destructor
426  , charac); // fields
427  }
428 
429  auto death_key_row = transverse_idx_.find(death_key); // Find the beginning of the row.
430  std::pair<typename Cam::iterator, bool> result_insert_cam;
431 
432  auto row_cell_it = death_key_row->second.row_->begin();
433 
434  while (row_cell_it != death_key_row->second.row_->end()) { // Traverse all cells in
435  // the row at index death_key.
436  Arith_element w = coeff_field_.times_minus(inv_x, row_cell_it->coefficient_);
437 
438  if (w != coeff_field_.additive_identity()) {
439  Column * curr_col = row_cell_it->self_col_;
440  ++row_cell_it;
441  // Disconnect the column from the rows in the CAM.
442  for (auto& col_cell : curr_col->col_) {
443  col_cell.base_hook_cam_h::unlink();
444  }
445 
446  // Remove the column from the CAM before modifying its value
447  cam_.erase(cam_.iterator_to(*curr_col));
448  // Proceed to the reduction of the column
449  plus_equal_column(*curr_col, a_ds, w);
450 
451  if (curr_col->col_.empty()) { // If the column is null
452  ds_repr_[curr_col->class_key_] = NULL;
453  column_pool_.destroy(curr_col); // delete curr_col;
454  } else {
455  // Find whether the column obtained is already in the CAM
456  result_insert_cam = cam_.insert(*curr_col);
457  if (result_insert_cam.second) { // If it was not in the CAM before: insertion has succeeded
458  for (auto& col_cell : curr_col->col_) {
459  // re-establish the row links
460  transverse_idx_[col_cell.key_].row_->push_front(col_cell);
461  }
462  } else { // There is already an identical column in the CAM:
463  // merge two disjoint sets.
464  dsets_.link(curr_col->class_key_,
465  result_insert_cam.first->class_key_);
466 
467  Simplex_key key_tmp = dsets_.find_set(curr_col->class_key_);
468  ds_repr_[key_tmp] = &(*(result_insert_cam.first));
469  result_insert_cam.first->class_key_ = key_tmp;
470  // intrusive containers don't own their elements, we have to release them manually
471  curr_col->col_.clear_and_dispose([&](Cell*p){cell_pool_.destroy(p);});
472  column_pool_.destroy(curr_col); // delete curr_col;
473  }
474  }
475  } else {
476  ++row_cell_it;
477  } // If w == 0, pass.
478  }
479 
480  // Because it is a killer simplex, set the data of sigma to null_key().
481  if (charac == coeff_field_.characteristic()) {
482  cpx_->assign_key(sigma, cpx_->null_key());
483  }
484  if (death_key_row->second.characteristics_ == charac) {
485  delete death_key_row->second.row_;
486  transverse_idx_.erase(death_key_row);
487  } else {
488  death_key_row->second.characteristics_ /= charac;
489  }
490  }
491 
492  /*
493  * Assign: target <- target + w * other.
494  */
495  void plus_equal_column(Column & target, A_ds_type const& other // value_type is pair<Simplex_key,Arith_element>
496  , Arith_element w) {
497  auto target_it = target.col_.begin();
498  auto other_it = other.begin();
499  while (target_it != target.col_.end() && other_it != other.end()) {
500  if (target_it->key_ < other_it->first) {
501  ++target_it;
502  } else {
503  if (target_it->key_ > other_it->first) {
504  Cell * cell_tmp = cell_pool_.construct(Cell(other_it->first // key
505  , coeff_field_.additive_identity(), &target));
506 
507  cell_tmp->coefficient_ = coeff_field_.plus_times_equal(cell_tmp->coefficient_, other_it->second, w);
508 
509  target.col_.insert(target_it, *cell_tmp);
510 
511  ++other_it;
512  } else { // it1->key == it2->key
513  // target_it->coefficient_ <- target_it->coefficient_ + other_it->second * w
514  target_it->coefficient_ = coeff_field_.plus_times_equal(target_it->coefficient_, other_it->second, w);
515  if (target_it->coefficient_ == coeff_field_.additive_identity()) {
516  auto tmp_it = target_it;
517  ++target_it;
518  ++other_it; // iterators remain valid
519  Cell * tmp_cell_ptr = &(*tmp_it);
520  target.col_.erase(tmp_it); // removed from column
521 
522  cell_pool_.destroy(tmp_cell_ptr); // delete from memory
523  } else {
524  ++target_it;
525  ++other_it;
526  }
527  }
528  }
529  }
530  while (other_it != other.end()) {
531  Cell * cell_tmp = cell_pool_.construct(Cell(other_it->first, coeff_field_.additive_identity(), &target));
532  cell_tmp->coefficient_ = coeff_field_.plus_times_equal(cell_tmp->coefficient_, other_it->second, w);
533  target.col_.insert(target.col_.end(), *cell_tmp);
534 
535  ++other_it;
536  }
537  }
538 
539  /*
540  * Compare two intervals by length.
541  */
542  struct cmp_intervals_by_length {
543  explicit cmp_intervals_by_length(FilteredComplex * sc)
544  : sc_(sc) {
545  }
546  bool operator()(const Persistent_interval & p1, const Persistent_interval & p2) {
547  return (sc_->filtration(get < 1 > (p1)) - sc_->filtration(get < 0 > (p1))
548  > sc_->filtration(get < 1 > (p2)) - sc_->filtration(get < 0 > (p2)));
549  }
550  FilteredComplex * sc_;
551  };
552 
553  public:
564  void output_diagram(std::ostream& ostream = std::cout) {
565  cmp_intervals_by_length cmp(cpx_);
566  std::sort(std::begin(persistent_pairs_), std::end(persistent_pairs_), cmp);
567  bool has_infinity = std::numeric_limits<Filtration_value>::has_infinity;
568  for (auto pair : persistent_pairs_) {
569  // Special case on windows, inf is "1.#INF" (cf. unitary tests and R package TDA)
570  if (has_infinity && cpx_->filtration(get<1>(pair)) == std::numeric_limits<Filtration_value>::infinity()) {
571  ostream << get<2>(pair) << " " << cpx_->dimension(get<0>(pair)) << " "
572  << cpx_->filtration(get<0>(pair)) << " inf " << std::endl;
573  } else {
574  ostream << get<2>(pair) << " " << cpx_->dimension(get<0>(pair)) << " "
575  << cpx_->filtration(get<0>(pair)) << " "
576  << cpx_->filtration(get<1>(pair)) << " " << std::endl;
577  }
578  }
579  }
580 
581  void write_output_diagram(std::string diagram_name) {
582  std::ofstream diagram_out(diagram_name.c_str());
583  cmp_intervals_by_length cmp(cpx_);
584  std::sort(std::begin(persistent_pairs_), std::end(persistent_pairs_), cmp);
585  bool has_infinity = std::numeric_limits<Filtration_value>::has_infinity;
586  for (auto pair : persistent_pairs_) {
587  // Special case on windows, inf is "1.#INF"
588  if (has_infinity && cpx_->filtration(get<1>(pair)) == std::numeric_limits<Filtration_value>::infinity()) {
589  diagram_out << cpx_->dimension(get<0>(pair)) << " "
590  << cpx_->filtration(get<0>(pair)) << " inf" << std::endl;
591  } else {
592  diagram_out << cpx_->dimension(get<0>(pair)) << " "
593  << cpx_->filtration(get<0>(pair)) << " "
594  << cpx_->filtration(get<1>(pair)) << std::endl;
595  }
596  }
597  }
598 
602  std::vector<int> betti_numbers() const {
603  // Don't allocate a vector of negative size for an empty complex
604  int siz = std::max(dim_max_, 0);
605  // Init Betti numbers vector with zeros until Simplicial complex dimension
606  std::vector<int> betti_numbers(siz);
607 
608  for (auto pair : persistent_pairs_) {
609  // Count never ended persistence intervals
610  if (cpx_->null_simplex() == get<1>(pair)) {
611  // Increment corresponding betti number
612  betti_numbers[cpx_->dimension(get<0>(pair))] += 1;
613  }
614  }
615  return betti_numbers;
616  }
617 
623  int betti_number(int dimension) const {
624  int betti_number = 0;
625 
626  for (auto pair : persistent_pairs_) {
627  // Count never ended persistence intervals
628  if (cpx_->null_simplex() == get<1>(pair)) {
629  if (cpx_->dimension(get<0>(pair)) == dimension) {
630  // Increment betti number found
631  ++betti_number;
632  }
633  }
634  }
635  return betti_number;
636  }
637 
643  std::vector<int> persistent_betti_numbers(Filtration_value from, Filtration_value to) const {
644  // Don't allocate a vector of negative size for an empty complex
645  int siz = std::max(dim_max_, 0);
646  // Init Betti numbers vector with zeros until Simplicial complex dimension
647  std::vector<int> betti_numbers(siz);
648  for (auto pair : persistent_pairs_) {
649  // Count persistence intervals that covers the given interval
650  // null_simplex test : if the function is called with to=+infinity, we still get something useful. And it will
651  // still work if we change the complex filtration function to reject null simplices.
652  if (cpx_->filtration(get<0>(pair)) <= from &&
653  (get<1>(pair) == cpx_->null_simplex() || cpx_->filtration(get<1>(pair)) > to)) {
654  // Increment corresponding betti number
655  betti_numbers[cpx_->dimension(get<0>(pair))] += 1;
656  }
657  }
658  return betti_numbers;
659  }
660 
667  int persistent_betti_number(int dimension, Filtration_value from, Filtration_value to) const {
668  int betti_number = 0;
669 
670  for (auto pair : persistent_pairs_) {
671  // Count persistence intervals that covers the given interval
672  // null_simplex test : if the function is called with to=+infinity, we still get something useful. And it will
673  // still work if we change the complex filtration function to reject null simplices.
674  if (cpx_->filtration(get<0>(pair)) <= from &&
675  (get<1>(pair) == cpx_->null_simplex() || cpx_->filtration(get<1>(pair)) > to)) {
676  if (cpx_->dimension(get<0>(pair)) == dimension) {
677  // Increment betti number found
678  ++betti_number;
679  }
680  }
681  }
682  return betti_number;
683  }
684 
688  const std::vector<Persistent_interval>& get_persistent_pairs() const {
689  return persistent_pairs_;
690  }
691 
696  std::vector< std::pair< Filtration_value , Filtration_value > >
697  intervals_in_dimension(int dimension) {
698  std::vector< std::pair< Filtration_value , Filtration_value > > result;
699  // auto && pair, to avoid unnecessary copying
700  for (auto && pair : persistent_pairs_) {
701  if (cpx_->dimension(get<0>(pair)) == dimension) {
702  result.emplace_back(cpx_->filtration(get<0>(pair)), cpx_->filtration(get<1>(pair)));
703  }
704  }
705  return result;
706  }
707 
708  private:
709  /*
710  * Structure representing a cocycle.
711  */
712  struct cocycle {
713  cocycle()
714  : row_(nullptr),
715  characteristics_() {
716  }
717  cocycle(Arith_element characteristics, Hcell * row)
718  : row_(row),
719  characteristics_(characteristics) {
720  }
721 
722  Hcell * row_; // points to the corresponding row in the CAM
723  Arith_element characteristics_; // product of field characteristics for which the cocycle exist
724  };
725 
726  public:
727  FilteredComplex * cpx_;
728  int dim_max_;
729  CoefficientField coeff_field_;
730  size_t num_simplices_;
731 
732  /* Disjoint sets data structure to link the model of FilteredComplex
733  * with the compressed annotation matrix.
734  * ds_rank_ is a property map Simplex_key -> int, ds_parent_ is a property map
735  * Simplex_key -> simplex_key_t */
736  std::vector<int> ds_rank_;
737  std::vector<Simplex_key> ds_parent_;
738  std::vector<Column *> ds_repr_;
739  boost::disjoint_sets<int *, Simplex_key *> dsets_;
740  /* The compressed annotation matrix fields.*/
741  Cam cam_;
742  /* Dictionary establishing the correspondance between the Simplex_key of
743  * the root vertex in the union-find ds and the Simplex_key of the vertex which
744  * created the connected component as a 0-dimension homology feature.*/
745  std::map<Simplex_key, Simplex_key> zero_cocycles_;
746  /* Key -> row. */
747  std::map<Simplex_key, cocycle> transverse_idx_;
748  /* Persistent intervals. */
749  std::vector<Persistent_interval> persistent_pairs_;
750  length_interval interval_length_policy;
751 
752  Simple_object_pool<Column> column_pool_;
753  Simple_object_pool<Cell> cell_pool_;
754 };
755 
756 } // namespace persistent_cohomology
757 
758 } // namespace Gudhi
759 
760 #endif // PERSISTENT_COHOMOLOGY_H_
FilteredComplex::Simplex_key Simplex_key
Data stored for each simplex.
Definition: Persistent_cohomology.h:57
void init_coefficients(int charac)
Initializes the coefficient field.
Definition: Persistent_cohomology.h:156
unspecified Element
Type of element of the field.
Definition: CoefficientField.h:19
unspecified Simplex_key
Data stored for each simplex.
Definition: FilteredComplex.h:91
int persistent_betti_number(int dimension, Filtration_value from, Filtration_value to) const
Returns the persistent Betti number of the dimension passed by parameter.
Definition: Persistent_cohomology.h:667
void assign_key(Simplex_handle sh, Simplex_key n)
Store a number for a simplex, which can later be retrieved with key(sh).
Computes the persistent cohomology of a filtered complex.
Definition: Persistent_cohomology.h:52
int betti_number(int dimension) const
Returns the Betti number of the dimension passed by parameter.
Definition: Persistent_cohomology.h:623
CoefficientField::Element Arith_element
Type of element of the field.
Definition: Persistent_cohomology.h:63
Definition: SimplicialComplexForAlpha.h:14
const std::vector< Persistent_interval > & get_persistent_pairs() const
Returns a list of persistence birth and death FilteredComplex::Simplex_handle pairs.
Definition: Persistent_cohomology.h:688
Simplex_key key(Simplex_handle sh)
Returns the number stored for a simplex by assign_key.
Persistent_cohomology(FilteredComplex &cpx, bool persistence_dim_max=false)
Initializes the Persistent_cohomology class.
Definition: Persistent_cohomology.h:95
std::vector< std::pair< Filtration_value, Filtration_value > > intervals_in_dimension(int dimension)
Returns persistence intervals for a given dimension.
Definition: Persistent_cohomology.h:697
void output_diagram(std::ostream &ostream=std::cout)
Output the persistence diagram in ostream.
Definition: Persistent_cohomology.h:564
Filtration_simplex_range filtration_simplex_range()
Returns a range over the simplices of the complex in the order of the filtration. ...
void init_coefficients(int charac_min, int charac_max)
Initializes the coefficient field for multi-field persistent homology.
Definition: Persistent_cohomology.h:160
std::tuple< Simplex_handle, Simplex_handle, Arith_element > Persistent_interval
Type for birth and death FilteredComplex::Simplex_handle. The Arith_element field is used for the mul...
Definition: Persistent_cohomology.h:66
size_t num_simplices()
Returns the number of simplices in the complex.
unspecified Simplex_handle
Handle to specify a simplex.
Definition: FilteredComplex.h:19
FilteredComplex::Simplex_handle Simplex_handle
Handle to specify a simplex.
Definition: Persistent_cohomology.h:59
std::vector< int > persistent_betti_numbers(Filtration_value from, Filtration_value to) const
Returns the persistent Betti numbers.
Definition: Persistent_cohomology.h:643
Concept describing the requirements for a class to represent a field of coefficients to compute persi...
Definition: CoefficientField.h:14
std::vector< int > betti_numbers() const
Returns Betti numbers.
Definition: Persistent_cohomology.h:602
FilteredComplex::Filtration_value Filtration_value
Type for the value of the filtration function.
Definition: Persistent_cohomology.h:61
void compute_persistent_cohomology(Filtration_value min_interval_length=0)
Compute the persistent homology of the filtered simplicial complex.
Definition: Persistent_cohomology.h:172
The concept FilteredComplex describes the requirements for a type to implement a filtered cell comple...
Definition: FilteredComplex.h:16
unspecified Filtration_value
Type for the value of the filtration function.
Definition: FilteredComplex.h:23
GUDHI  Version 3.1.1  - C++ library for Topological Data Analysis (TDA) and Higher Dimensional Geometry Understanding.  - Copyright : MIT Generated on Fri Feb 7 2020 16:35:36 for GUDHI by Doxygen 1.8.13