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
34namespace Gudhi {
35
36namespace persistent_cohomology {
37
50// TODO(CM): Memory allocation policy: classic, use a mempool, etc.
51template<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_.data(), ds_parent_.data()), // 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 thread_local std::vector<annotation_t> annotations_in_boundary;
292 annotations_in_boundary.clear();
293 int sign = 1 - 2 * (dim_sigma % 2); // \in {-1,1} provides the sign in the
294 // alternate sum in the boundary.
295 Simplex_key key;
296 Column * curr_col;
297
298 for (auto sh : cpx_->boundary_simplex_range(sigma)) {
299 key = cpx_->key(sh);
300 if (key != cpx_->null_key()) { // A simplex with null_key is a killer, and have null annotation
301 // Find its annotation vector
302 curr_col = ds_repr_[dsets_.find_set(key)];
303 if (curr_col != NULL) { // and insert it in annotations_in_boundary with multyiplicative factor "sign".
304 annotations_in_boundary.emplace_back(curr_col, sign);
305 }
306 }
307 sign = -sign;
308 }
309 // Place identical annotations consecutively so we can easily sum their multiplicities.
310 std::sort(annotations_in_boundary.begin(), annotations_in_boundary.end(),
311 [](annotation_t const& a, annotation_t const& b) { return a.first < b.first; });
312
313 // Sum the annotations with multiplicity, using a map<key,coeff>
314 // to represent a sparse vector.
315 std::pair<typename std::map<Simplex_key, Arith_element>::iterator, bool> result_insert_a_ds;
316
317 for (auto ann_it = annotations_in_boundary.begin(); ann_it != annotations_in_boundary.end(); ) {
318 Column* col = ann_it->first;
319 int mult = ann_it->second;
320 while (++ann_it != annotations_in_boundary.end() && ann_it->first == col) {
321 mult += ann_it->second;
322 }
323 // The following test is just a heuristic, it is not required, and it is fine that is misses p == 0.
324 if (mult != coeff_field_.additive_identity()) { // For all columns in the boundary,
325 for (auto cell_ref : col->col_) { // insert every cell in map_a_ds with multiplicity
326 Arith_element w_y = coeff_field_.times(cell_ref.coefficient_, mult); // coefficient * multiplicity
327
328 if (w_y != coeff_field_.additive_identity()) { // if != 0
329 result_insert_a_ds = map_a_ds.insert(std::pair<Simplex_key, Arith_element>(cell_ref.key_, w_y));
330 if (!(result_insert_a_ds.second)) { // if cell_ref.key_ already a Key in map_a_ds
331 result_insert_a_ds.first->second = coeff_field_.plus_equal(result_insert_a_ds.first->second, w_y);
332 if (result_insert_a_ds.first->second == coeff_field_.additive_identity()) {
333 map_a_ds.erase(result_insert_a_ds.first);
334 }
335 }
336 }
337 }
338 }
339 }
340 }
341
342 /*
343 * Update the cohomology groups under the insertion of a simplex.
344 */
345 void update_cohomology_groups(Simplex_handle sigma, int dim_sigma) {
346// Compute the annotation of the boundary of sigma:
347 std::map<Simplex_key, Arith_element> map_a_ds;
348 annotation_of_the_boundary(map_a_ds, sigma, dim_sigma);
349// Update the cohomology groups:
350 if (map_a_ds.empty()) { // sigma is a creator in all fields represented in coeff_field_
351 if (dim_sigma < dim_max_) {
352 create_cocycle(sigma, coeff_field_.multiplicative_identity(),
353 coeff_field_.characteristic());
354 }
355 } else { // sigma is a destructor in at least a field in coeff_field_
356 // Convert map_a_ds to a vector
357 A_ds_type a_ds; // admits reverse iterators
358 for (auto map_a_ds_ref : map_a_ds) {
359 a_ds.push_back(
360 std::pair<Simplex_key, Arith_element>(map_a_ds_ref.first,
361 map_a_ds_ref.second));
362 }
363
364 Arith_element inv_x, charac;
365 Arith_element prod = coeff_field_.characteristic(); // Product of characteristic of the fields
366 for (auto a_ds_rit = a_ds.rbegin();
367 (a_ds_rit != a_ds.rend())
368 && (prod != coeff_field_.multiplicative_identity()); ++a_ds_rit) {
369 std::tie(inv_x, charac) = coeff_field_.inverse(a_ds_rit->second, prod);
370
371 if (inv_x != coeff_field_.additive_identity()) {
372 destroy_cocycle(sigma, a_ds, a_ds_rit->first, inv_x, charac);
373 prod /= charac;
374 }
375 }
376 if (prod != coeff_field_.multiplicative_identity()
377 && dim_sigma < dim_max_) {
378 create_cocycle(sigma, coeff_field_.multiplicative_identity(prod), prod);
379 }
380 }
381 }
382
383 /* \brief Create a new cocycle class.
384 *
385 * The class is created by the insertion of the simplex sigma.
386 * The methods adds a cocycle, representing the new cocycle class,
387 * to the matrix representing the cohomology groups.
388 * The new cocycle has value 0 on every simplex except on sigma
389 * where it worths 1.*/
390 void create_cocycle(Simplex_handle sigma, Arith_element x,
391 Arith_element charac) {
392 Simplex_key key = cpx_->key(sigma);
393 // Create a column containing only one cell,
394 Column * new_col = column_pool_.construct(key);
395 Cell * new_cell = cell_pool_.construct(key, x, new_col);
396 new_col->col_.push_back(*new_cell);
397 // and insert it in the matrix, in constant time thanks to the hint cam_.end().
398 // Indeed *new_col has the biggest lexicographic value because key is the
399 // biggest key used so far.
400 cam_.insert(cam_.end(), *new_col);
401 // Update the disjoint sets data structure.
402 Hcell * new_hcell = new Hcell;
403 new_hcell->push_back(*new_cell);
404 transverse_idx_[key] = cocycle(charac, new_hcell); // insert the new row
405 ds_repr_[key] = new_col;
406 }
407
408 /* \brief Destroy a cocycle class.
409 *
410 * The cocycle class is destroyed by the insertion of sigma.
411 * The methods proceeds to a reduction of the matrix representing
412 * the cohomology groups using Gauss pivoting. The reduction zeros-out
413 * the row containing the cell with highest key in
414 * a_ds, the annotation of the boundary of simplex sigma. This key
415 * is "death_key".*/
416 void destroy_cocycle(Simplex_handle sigma, A_ds_type const& a_ds,
417 Simplex_key death_key, Arith_element inv_x,
418 Arith_element charac) {
419 // Create a finite persistent interval for which the interval exists
420 if (interval_length_policy(cpx_->simplex(death_key), sigma)) {
421 persistent_pairs_.emplace_back(cpx_->simplex(death_key) // creator
422 , sigma // destructor
423 , charac); // fields
424 }
425
426 auto death_key_row = transverse_idx_.find(death_key); // Find the beginning of the row.
427 std::pair<typename Cam::iterator, bool> result_insert_cam;
428
429 auto row_cell_it = death_key_row->second.row_->begin();
430
431 while (row_cell_it != death_key_row->second.row_->end()) { // Traverse all cells in
432 // the row at index death_key.
433 Arith_element w = coeff_field_.times_minus(inv_x, row_cell_it->coefficient_);
434
435 if (w != coeff_field_.additive_identity()) {
436 Column * curr_col = row_cell_it->self_col_;
437 ++row_cell_it;
438 // Disconnect the column from the rows in the CAM.
439 for (auto& col_cell : curr_col->col_) {
440 col_cell.base_hook_cam_h::unlink();
441 }
442
443 // Remove the column from the CAM before modifying its value
444 cam_.erase(cam_.iterator_to(*curr_col));
445 // Proceed to the reduction of the column
446 plus_equal_column(*curr_col, a_ds, w);
447
448 if (curr_col->col_.empty()) { // If the column is null
449 ds_repr_[curr_col->class_key_] = NULL;
450 column_pool_.destroy(curr_col); // delete curr_col;
451 } else {
452 // Find whether the column obtained is already in the CAM
453 result_insert_cam = cam_.insert(*curr_col);
454 if (result_insert_cam.second) { // If it was not in the CAM before: insertion has succeeded
455 for (auto& col_cell : curr_col->col_) {
456 // re-establish the row links
457 transverse_idx_[col_cell.key_].row_->push_front(col_cell);
458 }
459 } else { // There is already an identical column in the CAM:
460 // merge two disjoint sets.
461 dsets_.link(curr_col->class_key_,
462 result_insert_cam.first->class_key_);
463
464 Simplex_key key_tmp = dsets_.find_set(curr_col->class_key_);
465 ds_repr_[key_tmp] = &(*(result_insert_cam.first));
466 result_insert_cam.first->class_key_ = key_tmp;
467 // intrusive containers don't own their elements, we have to release them manually
468 curr_col->col_.clear_and_dispose([&](Cell*p){cell_pool_.destroy(p);});
469 column_pool_.destroy(curr_col); // delete curr_col;
470 }
471 }
472 } else {
473 ++row_cell_it;
474 } // If w == 0, pass.
475 }
476
477 // Because it is a killer simplex, set the data of sigma to null_key().
478 if (charac == coeff_field_.characteristic()) {
479 cpx_->assign_key(sigma, cpx_->null_key());
480 }
481 if (death_key_row->second.characteristics_ == charac) {
482 delete death_key_row->second.row_;
483 transverse_idx_.erase(death_key_row);
484 } else {
485 death_key_row->second.characteristics_ /= charac;
486 }
487 }
488
489 /*
490 * Assign: target <- target + w * other.
491 */
492 void plus_equal_column(Column & target, A_ds_type const& other // value_type is pair<Simplex_key,Arith_element>
493 , Arith_element w) {
494 auto target_it = target.col_.begin();
495 auto other_it = other.begin();
496 while (target_it != target.col_.end() && other_it != other.end()) {
497 if (target_it->key_ < other_it->first) {
498 ++target_it;
499 } else {
500 if (target_it->key_ > other_it->first) {
501 Cell * cell_tmp = cell_pool_.construct(Cell(other_it->first // key
502 , coeff_field_.additive_identity(), &target));
503
504 cell_tmp->coefficient_ = coeff_field_.plus_times_equal(cell_tmp->coefficient_, other_it->second, w);
505
506 target.col_.insert(target_it, *cell_tmp);
507
508 ++other_it;
509 } else { // it1->key == it2->key
510 // target_it->coefficient_ <- target_it->coefficient_ + other_it->second * w
511 target_it->coefficient_ = coeff_field_.plus_times_equal(target_it->coefficient_, other_it->second, w);
512 if (target_it->coefficient_ == coeff_field_.additive_identity()) {
513 auto tmp_it = target_it;
514 ++target_it;
515 ++other_it; // iterators remain valid
516 Cell * tmp_cell_ptr = &(*tmp_it);
517 target.col_.erase(tmp_it); // removed from column
518
519 cell_pool_.destroy(tmp_cell_ptr); // delete from memory
520 } else {
521 ++target_it;
522 ++other_it;
523 }
524 }
525 }
526 }
527 while (other_it != other.end()) {
528 Cell * cell_tmp = cell_pool_.construct(Cell(other_it->first, coeff_field_.additive_identity(), &target));
529 cell_tmp->coefficient_ = coeff_field_.plus_times_equal(cell_tmp->coefficient_, other_it->second, w);
530 target.col_.insert(target.col_.end(), *cell_tmp);
531
532 ++other_it;
533 }
534 }
535
536 /*
537 * Compare two intervals by length.
538 */
539 struct cmp_intervals_by_length {
540 explicit cmp_intervals_by_length(FilteredComplex * sc)
541 : sc_(sc) {
542 }
543 bool operator()(const Persistent_interval & p1, const Persistent_interval & p2) {
544 return (sc_->filtration(get < 1 > (p1)) - sc_->filtration(get < 0 > (p1))
545 > sc_->filtration(get < 1 > (p2)) - sc_->filtration(get < 0 > (p2)));
546 }
547 FilteredComplex * sc_;
548 };
549
550 public:
561 void output_diagram(std::ostream& ostream = std::cout) {
562 cmp_intervals_by_length cmp(cpx_);
563 std::sort(std::begin(persistent_pairs_), std::end(persistent_pairs_), cmp);
564 for (auto pair : persistent_pairs_) {
565 ostream << get<2>(pair) << " " << cpx_->dimension(get<0>(pair)) << " "
566 << cpx_->filtration(get<0>(pair)) << " "
567 << cpx_->filtration(get<1>(pair)) << " " << std::endl;
568 }
569 }
570
571 void write_output_diagram(std::string diagram_name) {
572 std::ofstream diagram_out(diagram_name.c_str());
573 diagram_out.exceptions(diagram_out.failbit);
574 cmp_intervals_by_length cmp(cpx_);
575 std::sort(std::begin(persistent_pairs_), std::end(persistent_pairs_), cmp);
576 for (auto pair : persistent_pairs_) {
577 diagram_out << cpx_->dimension(get<0>(pair)) << " "
578 << cpx_->filtration(get<0>(pair)) << " "
579 << cpx_->filtration(get<1>(pair)) << std::endl;
580 }
581 }
582
586 std::vector<int> betti_numbers() const {
587 // Don't allocate a vector of negative size for an empty complex
588 int siz = std::max(dim_max_, 0);
589 // Init Betti numbers vector with zeros until Simplicial complex dimension
590 std::vector<int> betti_numbers(siz);
591
592 for (auto pair : persistent_pairs_) {
593 // Count never ended persistence intervals
594 if (cpx_->null_simplex() == get<1>(pair)) {
595 // Increment corresponding betti number
596 betti_numbers[cpx_->dimension(get<0>(pair))] += 1;
597 }
598 }
599 return betti_numbers;
600 }
601
607 int betti_number(int dimension) const {
608 int betti_number = 0;
609
610 for (auto pair : persistent_pairs_) {
611 // Count never ended persistence intervals
612 if (cpx_->null_simplex() == get<1>(pair)) {
613 if (cpx_->dimension(get<0>(pair)) == dimension) {
614 // Increment betti number found
615 ++betti_number;
616 }
617 }
618 }
619 return betti_number;
620 }
621
628 // Don't allocate a vector of negative size for an empty complex
629 int siz = std::max(dim_max_, 0);
630 // Init Betti numbers vector with zeros until Simplicial complex dimension
631 std::vector<int> betti_numbers(siz);
632 for (auto pair : persistent_pairs_) {
633 // Count persistence intervals that covers the given interval
634 // null_simplex test : if the function is called with to=+infinity, we still get something useful. And it will
635 // still work if we change the complex filtration function to reject null simplices.
636 if (cpx_->filtration(get<0>(pair)) <= from &&
637 (get<1>(pair) == cpx_->null_simplex() || cpx_->filtration(get<1>(pair)) > to)) {
638 // Increment corresponding betti number
639 betti_numbers[cpx_->dimension(get<0>(pair))] += 1;
640 }
641 }
642 return betti_numbers;
643 }
644
651 int persistent_betti_number(int dimension, Filtration_value from, Filtration_value to) const {
652 int betti_number = 0;
653
654 for (auto pair : persistent_pairs_) {
655 // Count persistence intervals that covers the given interval
656 // null_simplex test : if the function is called with to=+infinity, we still get something useful. And it will
657 // still work if we change the complex filtration function to reject null simplices.
658 if (cpx_->filtration(get<0>(pair)) <= from &&
659 (get<1>(pair) == cpx_->null_simplex() || cpx_->filtration(get<1>(pair)) > to)) {
660 if (cpx_->dimension(get<0>(pair)) == dimension) {
661 // Increment betti number found
662 ++betti_number;
663 }
664 }
665 }
666 return betti_number;
667 }
668
672 const std::vector<Persistent_interval>& get_persistent_pairs() const {
673 return persistent_pairs_;
674 }
675
680 std::vector< std::pair< Filtration_value , Filtration_value > >
681 intervals_in_dimension(int dimension) {
682 std::vector< std::pair< Filtration_value , Filtration_value > > result;
683 // auto && pair, to avoid unnecessary copying
684 for (auto && pair : persistent_pairs_) {
685 if (cpx_->dimension(get<0>(pair)) == dimension) {
686 result.emplace_back(cpx_->filtration(get<0>(pair)), cpx_->filtration(get<1>(pair)));
687 }
688 }
689 return result;
690 }
691
692 private:
693 /*
694 * Structure representing a cocycle.
695 */
696 struct cocycle {
697 cocycle()
698 : row_(nullptr),
699 characteristics_() {
700 }
701 cocycle(Arith_element characteristics, Hcell * row)
702 : row_(row),
703 characteristics_(characteristics) {
704 }
705
706 Hcell * row_; // points to the corresponding row in the CAM
707 Arith_element characteristics_; // product of field characteristics for which the cocycle exist
708 };
709
710 public:
711 FilteredComplex * cpx_;
712 int dim_max_;
713 CoefficientField coeff_field_;
714 size_t num_simplices_;
715
716 /* Disjoint sets data structure to link the model of FilteredComplex
717 * with the compressed annotation matrix.
718 * ds_rank_ is a property map Simplex_key -> int, ds_parent_ is a property map
719 * Simplex_key -> simplex_key_t */
720 std::vector<int> ds_rank_;
721 std::vector<Simplex_key> ds_parent_;
722 std::vector<Column *> ds_repr_;
723 boost::disjoint_sets<int *, Simplex_key *> dsets_;
724 /* The compressed annotation matrix fields.*/
725 Cam cam_;
726 /* Dictionary establishing the correspondance between the Simplex_key of
727 * the root vertex in the union-find ds and the Simplex_key of the vertex which
728 * created the connected component as a 0-dimension homology feature.*/
729 std::map<Simplex_key, Simplex_key> zero_cocycles_;
730 /* Key -> row. */
731 std::map<Simplex_key, cocycle> transverse_idx_;
732 /* Persistent intervals. */
733 std::vector<Persistent_interval> persistent_pairs_;
734 length_interval interval_length_policy;
735
736 Simple_object_pool<Column> column_pool_;
737 Simple_object_pool<Cell> cell_pool_;
738};
739
740} // namespace persistent_cohomology
741
742} // namespace Gudhi
743
744#endif // PERSISTENT_COHOMOLOGY_H_
Computes the persistent cohomology of a filtered complex.
Definition: Persistent_cohomology.h:52
std::vector< int > persistent_betti_numbers(Filtration_value from, Filtration_value to) const
Returns the persistent Betti numbers.
Definition: Persistent_cohomology.h:627
std::vector< std::pair< Filtration_value, Filtration_value > > intervals_in_dimension(int dimension)
Returns persistence intervals for a given dimension.
Definition: Persistent_cohomology.h:681
std::vector< int > betti_numbers() const
Returns Betti numbers.
Definition: Persistent_cohomology.h:586
void output_diagram(std::ostream &ostream=std::cout)
Output the persistence diagram in ostream.
Definition: Persistent_cohomology.h:561
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
FilteredComplex::Simplex_handle Simplex_handle
Handle to specify a simplex.
Definition: Persistent_cohomology.h:59
Persistent_cohomology(FilteredComplex &cpx, bool persistence_dim_max=false)
Initializes the Persistent_cohomology class.
Definition: Persistent_cohomology.h:95
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:651
int betti_number(int dimension) const
Returns the Betti number of the dimension passed by parameter.
Definition: Persistent_cohomology.h:607
FilteredComplex::Simplex_key Simplex_key
Data stored for each simplex.
Definition: Persistent_cohomology.h:57
void compute_persistent_cohomology(Filtration_value min_interval_length=0)
Compute the persistent homology of the filtered simplicial complex.
Definition: Persistent_cohomology.h:172
void init_coefficients(int charac)
Initializes the coefficient field.
Definition: Persistent_cohomology.h:156
void init_coefficients(int charac_min, int charac_max)
Initializes the coefficient field for multi-field persistent homology.
Definition: Persistent_cohomology.h:160
CoefficientField::Element Arith_element
Type of element of the field.
Definition: Persistent_cohomology.h:63
FilteredComplex::Filtration_value Filtration_value
Type for the value of the filtration function.
Definition: Persistent_cohomology.h:61
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:672
Concept describing the requirements for a class to represent a field of coefficients to compute persi...
Definition: CoefficientField.h:14
Element additive_identity()
Element multiplicative_identity()
unspecified Element
Type of element of the field.
Definition: CoefficientField.h:19
Element characteristic()
void plus_equal(Element x, Element y)
The concept FilteredComplex describes the requirements for a type to implement a filtered cell comple...
Definition: FilteredComplex.h:17
unspecified Simplex_key
Data stored for each simplex.
Definition: FilteredComplex.h:91
Filtration_value filtration(Simplex_handle sh)
Returns the filtration value of a simplex.
void assign_key(Simplex_handle sh, Simplex_key n)
Store a number for a simplex, which can later be retrieved with key(sh).
Simplex_handle null_simplex()
Returns a Simplex_handle that is different from all simplex handles of the simplices.
Filtration_simplex_range filtration_simplex_range()
Returns a range over the simplices of the complex in the order of the filtration.
unspecified Simplex_handle
Handle to specify a simplex.
Definition: FilteredComplex.h:19
Simplex_key null_key()
Returns a constant dummy number that is either negative, or at least as large as num_simplices()....
Simplex_key key(Simplex_handle sh)
Returns the number stored for a simplex by assign_key.
unspecified Filtration_value
Type for the value of the filtration function.
Definition: FilteredComplex.h:23
int dimension(Simplex_handle sh)
Returns the dimension of a simplex.
Simplex_handle simplex(size_t idx)
Returns the simplex that has index idx in the filtration.
Boundary_simplex_range boundary_simplex_range(Simplex_handle sh)
Returns a range giving access to all simplices of the boundary of a simplex, i.e. the set of codimens...
size_t num_simplices()
Returns the number of simplices in the complex.
Value type for a filtration function on a cell complex.
Definition: FiltrationValue.h:20