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