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