11 #ifndef PERSISTENT_COHOMOLOGY_H_
12 #define PERSISTENT_COHOMOLOGY_H_
14 #include <gudhi/Persistent_cohomology/Persistent_cohomology_column.h>
15 #include <gudhi/Persistent_cohomology/Field_Zp.h>
16 #include <gudhi/Simple_object_pool.h>
18 #include <boost/intrusive/set.hpp>
19 #include <boost/pending/disjoint_sets.hpp>
20 #include <boost/intrusive/list.hpp>
36 namespace persistent_cohomology {
51 template<
class FilteredComplex,
class CoefficientField>
71 typedef Persistent_cohomology_column<Simplex_key, Arith_element> Column;
73 typedef typename Column::Cell Cell;
75 typedef boost::intrusive::list<Cell,
76 boost::intrusive::constant_time_size<false>,
77 boost::intrusive::base_hook<base_hook_cam_h> > Hcell;
79 typedef boost::intrusive::set<Column,
80 boost::intrusive::constant_time_size<false> > Cam;
82 typedef std::vector<std::pair<Simplex_key, Arith_element> > A_ds_type;
97 dim_max_(cpx.dimension()),
99 num_simplices_(cpx_->num_simplices()),
100 ds_rank_(num_simplices_),
101 ds_parent_(num_simplices_),
102 ds_repr_(num_simplices_, NULL),
103 dsets_(&ds_rank_[0], &ds_parent_[0]),
108 interval_length_policy(&cpx, 0),
111 if (cpx_->
num_simplices() > std::numeric_limits<Simplex_key>::max()) {
113 throw std::out_of_range(
"The number of simplices is more than Simplex_key type numeric limit.");
119 dsets_.make_set(cpx_->
key(sh));
121 if (persistence_dim_max) {
128 for (
auto & transverse_ref : transverse_idx_) {
130 transverse_ref.second.row_->clear_and_dispose([&](Cell*p){p->~Cell();});
131 delete transverse_ref.second.row_;
136 struct length_interval {
139 min_length_(min_length) {
143 return cpx_->filtration(sh2) - cpx_->filtration(sh1) > min_length_;
147 min_length_ = new_length;
157 coeff_field_.init(charac);
161 coeff_field_.init(charac_min, charac_max);
173 interval_length_policy.set_length(min_interval_length);
177 switch (dim_simplex) {
181 update_cohomology_groups_edge(sh);
184 update_cohomology_groups(sh, dim_simplex);
190 for (
auto v_sh : cpx_->skeleton_simplex_range(0)) {
191 key = cpx_->
key(v_sh);
193 if (ds_parent_[key] == key
194 && zero_cocycles_.find(key) == zero_cocycles_.end()) {
195 persistent_pairs_.emplace_back(
199 for (
auto zero_idx : zero_cocycles_) {
200 persistent_pairs_.emplace_back(
204 for (
auto cocycle : transverse_idx_) {
205 persistent_pairs_.emplace_back(
217 boost::tie(u, v) = cpx_->endpoints(sigma);
227 auto map_it_u = zero_cocycles_.find(ku);
229 if (map_it_u == zero_cocycles_.end()) {
232 idx_coc_u = map_it_u->second;
235 auto map_it_v = zero_cocycles_.find(kv);
237 if (map_it_v == zero_cocycles_.end()) {
240 idx_coc_v = map_it_v->second;
245 if (interval_length_policy(cpx_->
simplex(idx_coc_v), sigma)) {
246 persistent_pairs_.emplace_back(
250 if (kv != idx_coc_v) {
251 zero_cocycles_.erase(map_it_v);
253 if (kv == dsets_.find_set(kv)) {
254 if (ku != idx_coc_u) {
255 zero_cocycles_.erase(map_it_u);
257 zero_cocycles_[kv] = idx_coc_u;
260 if (interval_length_policy(cpx_->
simplex(idx_coc_u), sigma)) {
261 persistent_pairs_.emplace_back(
265 if (ku != idx_coc_u) {
266 zero_cocycles_.erase(map_it_u);
268 if (ku == dsets_.find_set(ku)) {
269 if (kv != idx_coc_v) {
270 zero_cocycles_.erase(map_it_v);
272 zero_cocycles_[ku] = idx_coc_v;
276 }
else if (dim_max_ > 1) {
284 void annotation_of_the_boundary(
285 std::map<Simplex_key, Arith_element> & map_a_ds,
Simplex_handle sigma,
290 typedef std::pair<Column *, int> annotation_t;
291 #ifdef GUDHI_CAN_USE_CXX11_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);
305 curr_col = ds_repr_[dsets_.find_set(key)];
306 if (curr_col != NULL) {
307 annotations_in_boundary.emplace_back(curr_col, sign);
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; });
318 std::pair<typename std::map<Simplex_key, Arith_element>::iterator,
bool> result_insert_a_ds;
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;
328 for (
auto cell_ref : col->col_) {
329 Arith_element w_y = coeff_field_.times(cell_ref.coefficient_, mult);
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)) {
334 result_insert_a_ds.first->second = coeff_field_.
plus_equal(result_insert_a_ds.first->second, w_y);
336 map_a_ds.erase(result_insert_a_ds.first);
348 void update_cohomology_groups(
Simplex_handle sigma,
int dim_sigma) {
350 std::map<Simplex_key, Arith_element> map_a_ds;
351 annotation_of_the_boundary(map_a_ds, sigma, dim_sigma);
353 if (map_a_ds.empty()) {
354 if (dim_sigma < dim_max_) {
361 for (
auto map_a_ds_ref : map_a_ds) {
363 std::pair<Simplex_key, Arith_element>(map_a_ds_ref.first,
364 map_a_ds_ref.second));
369 for (
auto a_ds_rit = a_ds.rbegin();
370 (a_ds_rit != a_ds.rend())
372 std::tie(inv_x, charac) = coeff_field_.inverse(a_ds_rit->second, prod);
375 destroy_cocycle(sigma, a_ds, a_ds_rit->first, inv_x, charac);
380 && dim_sigma < dim_max_) {
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);
403 cam_.insert(cam_.end(), *new_col);
405 Hcell * new_hcell =
new Hcell;
406 new_hcell->push_back(*new_cell);
407 transverse_idx_[key] = cocycle(charac, new_hcell);
408 ds_repr_[key] = new_col;
423 if (interval_length_policy(cpx_->
simplex(death_key), sigma)) {
424 persistent_pairs_.emplace_back(cpx_->
simplex(death_key)
429 auto death_key_row = transverse_idx_.find(death_key);
430 std::pair<typename Cam::iterator, bool> result_insert_cam;
432 auto row_cell_it = death_key_row->second.row_->begin();
434 while (row_cell_it != death_key_row->second.row_->end()) {
436 Arith_element w = coeff_field_.times_minus(inv_x, row_cell_it->coefficient_);
439 Column * curr_col = row_cell_it->self_col_;
442 for (
auto& col_cell : curr_col->col_) {
443 col_cell.base_hook_cam_h::unlink();
447 cam_.erase(cam_.iterator_to(*curr_col));
449 plus_equal_column(*curr_col, a_ds, w);
451 if (curr_col->col_.empty()) {
452 ds_repr_[curr_col->class_key_] = NULL;
453 column_pool_.destroy(curr_col);
456 result_insert_cam = cam_.insert(*curr_col);
457 if (result_insert_cam.second) {
458 for (
auto& col_cell : curr_col->col_) {
460 transverse_idx_[col_cell.key_].row_->push_front(col_cell);
464 dsets_.link(curr_col->class_key_,
465 result_insert_cam.first->class_key_);
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;
471 curr_col->col_.clear_and_dispose([&](Cell*p){cell_pool_.destroy(p);});
472 column_pool_.destroy(curr_col);
484 if (death_key_row->second.characteristics_ == charac) {
485 delete death_key_row->second.row_;
486 transverse_idx_.erase(death_key_row);
488 death_key_row->second.characteristics_ /= charac;
495 void plus_equal_column(Column & target, A_ds_type
const& other
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) {
503 if (target_it->key_ > other_it->first) {
504 Cell * cell_tmp = cell_pool_.construct(Cell(other_it->first
507 cell_tmp->coefficient_ = coeff_field_.plus_times_equal(cell_tmp->coefficient_, other_it->second, w);
509 target.col_.insert(target_it, *cell_tmp);
514 target_it->coefficient_ = coeff_field_.plus_times_equal(target_it->coefficient_, other_it->second, w);
516 auto tmp_it = target_it;
519 Cell * tmp_cell_ptr = &(*tmp_it);
520 target.col_.erase(tmp_it);
522 cell_pool_.destroy(tmp_cell_ptr);
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);
542 struct cmp_intervals_by_length {
547 return (sc_->filtration(get < 1 > (p1)) - sc_->filtration(get < 0 > (p1))
548 > sc_->filtration(get < 1 > (p2)) - sc_->filtration(get < 0 > (p2)));
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_) {
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;
574 ostream << get<2>(pair) <<
" " << cpx_->
dimension(get<0>(pair)) <<
" "
576 << cpx_->
filtration(get<1>(pair)) <<
" " << std::endl;
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_) {
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;
592 diagram_out << cpx_->
dimension(get<0>(pair)) <<
" "
594 << cpx_->
filtration(get<1>(pair)) << std::endl;
604 int siz = std::max(dim_max_, 0);
608 for (
auto pair : persistent_pairs_) {
626 for (
auto pair : persistent_pairs_) {
629 if (cpx_->
dimension(get<0>(pair)) == dimension) {
645 int siz = std::max(dim_max_, 0);
648 for (
auto pair : persistent_pairs_) {
670 for (
auto pair : persistent_pairs_) {
676 if (cpx_->
dimension(get<0>(pair)) == dimension) {
689 return persistent_pairs_;
696 std::vector< std::pair< Filtration_value , Filtration_value > >
698 std::vector< std::pair< Filtration_value , Filtration_value > > result;
700 for (
auto && pair : persistent_pairs_) {
701 if (cpx_->
dimension(get<0>(pair)) == dimension) {
719 characteristics_(characteristics) {
730 size_t num_simplices_;
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_;
745 std::map<Simplex_key, Simplex_key> zero_cocycles_;
747 std::map<Simplex_key, cocycle> transverse_idx_;
749 std::vector<Persistent_interval> persistent_pairs_;
750 length_interval interval_length_policy;
752 Simple_object_pool<Column> column_pool_;
753 Simple_object_pool<Cell> cell_pool_;
760 #endif // PERSISTENT_COHOMOLOGY_H_